The Utility of Mixed-Effect Models in the Evaluation of Complex Genomic Traits In-Vitro
[Methods Supplement and R Code]
Nathan Alade
January, 2023
Introduction
Abstract
The use of human liver microsomes as a model system to evaluate the impact of complex genomic traits (i.e., linkage-disequilibrium patterns, coding, and non-coding variation, etc.) on drug metabolism remains challenging. To overcome these challenges, we propose the use of non-linear mixed effect models (NLME) as an unconventional approach to evaluate the impact of complex genomic traits on the metabolism of xenobiotics in vitro. In this in-silico study, we present a practical use case for such an approach using previously published in-vitro CYP2D6 data and explore the impact of sparse sampling, and experimental error on known kinetic parameter of CYP2D6 mediated formation of 4-hydroxy-atomoxetine in human liver microsomes.
Objective
This supplement provides the R code used to conduct the
in-silico portion of the study, as well as useful functions for
evaluating/visualizing results.
R Libraries and Packages
# R Packages -----
library(nlme) # Non-Linear Mixed Effect Modeling
library(tidyverse) # Data Manipulation and Visualization
library(ggeasy) # Wrapper for Plot Manipulation
library(ggpubr) # Additional Tools for Plot Configuration
library(scales) # Additional Tools for Adjusting Plot Scales
library(cowplot) # Data Visualization Themes and Tools
library(kableExtra) # Additional Tools for Creating Custom Tables
# Package References -------
# knitr::write_bib(c(.packages()), "References/packages.bib")Citation information for R packages1–16 used for data analysis can be found in the References section.
Importing Data
CYP2D6 mean estimates and frequencies (Dinh et. al., 2016) were imported and summarized by genotype (referred to as diplotype in R script). Additionally, an Eta table and parameter estimate table were also generated to be reference later on in the data analysis.
# Atomoxetine Dataset (Dinh et al., 2016)
atx.data <- readxl::read_xlsx("Input Data/Atomoxetine Data - 2016.xlsx")
# Summary of Atomoxetine Data by genotype ----
diplotype.data <- atx.data %>%
group_by(Diplotype) %>%
# Summarize mean estimates of the data
summarise(Vmax = mean(Vmax),
Km = mean(Km),
Score = unique(`Activity Score`),
n = n()) %>%
ungroup() %>%
mutate(Freq = n/sum(n),
Vmax.Eta = ifelse(Diplotype != "1/1",
round(Vmax-Vmax[Diplotype == "1/1"],
digits = 2),
round(Vmax, digits = 2)),
Km.Eta = ifelse(Diplotype != "1/1",
round(Km-Km[Diplotype == "1/1"],
digits = 2),
round(Km, digits = 2))) %>%
select(-n)
# Eta Table (Deviations from the reference group) -----
actual.eta <- diplotype.data %>%
select(Diplotype, Vmax.Eta, Km.Eta) %>%
rename("Vmax" = "Vmax.Eta",
"Km" = "Km.Eta") %>%
pivot_longer(Vmax:Km, names_to = "Parameter",
values_to = "Actual") %>%
arrange(desc(Parameter))
# mean estimates of Vmax and Km for each genotype -----
actual.est <- diplotype.data %>%
select(Diplotype, Vmax, Km) %>%
pivot_longer(Vmax:Km, names_to = "Parameter",
values_to = "Actual") %>%
mutate(Actual = round(Actual, digits = 2)) %>%
arrange(desc(Parameter))Virtual Population
Population Simulation Function
A custom simulation function population.sim() was
scripted to generate a virtual population with a pre-defined
distribution of CYP2D6 genotypes, and to simulate
Michaelis-Menten kinetics for each individual. The inputs for the
function the population.sim() include the following:
n - number of subjects per genotype group (default = 10).
CV%_D - Average standard deviation for each genotype group as a percentage of the true estimate for Vmax and Km (between subject variability for each genotype group) (default = 25).
CV%_V - Average coefficient of variation across all points for simulated Michaelis-Menten kinetics for each individual (residual error).
Start - Starting concentration from Michaelis-Menten simulation (default = 0).
End - Ending concentration for Michaelis menten Kinetics (default = 100 uM).
By - Simulated concentration increment (default = 0.5).
pop_freq - TRUE or FALSE argument indicating whether to use true population frequency (default is FALSE).
seed - Random number generator seed (for reproducibility, default = 23457).
data - Dataframe to be used as reference for population averages. Must include genotype, Vmax, Km, and Frequency columns (default =
diplotype.data).
Within-group between subject variability was set at 25% CV for each genotype (one standard deviation = 0.25 x mean of true value). Below is a simulated population containing 9000 individuals (1000 per diplotype groups) with mean and standard deviations for each group parameter.
population.sim <- function(
# Pre-define number of subjects per diplotype group
n = 10,
`CV%_D` = 25,
`CV%_V` = 0,
# Pre-define Michaelis Menten Kinetic Settings
Start = 0,
End = 100,
By = 0.5,
# Turn off or on use of true population frequency
## If set to true "n" will equal total population number
pop_freq = FALSE,
# Pre-define reference data, and set seed
data = diplotype.data,
seed = 23456){
# ======== Create Data Containers ========
datasets <- list()
pop.data <- tibble()
# ======== Create Datasets by Diplotype ========
for (i in seq_along(data[[1]])){
# Setting Population Specific Frequency
n_pop <- if_else(pop_freq == FALSE, n,
ifelse(pop_freq == TRUE,
round(n*data[[i, c("Freq")]]), 0))
# Specify Diplotype for current loop
Diplotype = rep(data[[i,1]], n_pop)
# Random assignment of Vmax values N(0,sigma)
set.seed(seed)
Vmax_eta <- rnorm(n_pop, sd = `CV%_D`)
Vmax = data[[i,2]] + (Vmax_eta/100)*data[[i,2]]
# Random assignment of Km values N(0,sigma)
set.seed(seed)
Km_eta <- rnorm(n_pop, sd = `CV%_D`)
Km = data[[i,3]] + (Km_eta/100)*data[[i,3]]
# Store each diplotype group in a list container
temp <- tibble(Diplotype, Vmax, Km)
datasets[[i]] <- temp
}
# ======== Combine Diplotype Datasets ========
for (i in seq_along(data[[1]])){
pop.data <- rbind(pop.data, datasets[[i]])
}
# ======== Simulate Michaelis Menten ========
set.seed(seed)
pop.data <- pop.data %>%
mutate(ID = factor(seq_along(Diplotype),
levels = seq_along(Diplotype))) %>%
relocate(ID) %>%
expand_grid(S = seq(Start, End, By)) %>%
mutate(V = round((Vmax*S)/(Km+S), digits = 2))%>%
group_by(ID) %>%
# Additon of Residual Error ------------
mutate(resid_pct = round(rnorm(n = length(S),sd = `CV%_V`),
digits = 3),
V_adj = V+(V*resid_pct/100)) %>%
ungroup()
return(pop.data)
}Distribution of Virtual Population Parameters
dist.plot <- ggplot(
data = population.sim(n=1000,`CV%_D` = 25, Start = 0.5,
By = 5, pop_freq = F) %>%
group_by(ID) %>%
filter(S == min(S)),
aes(x = Vmax))+
geom_histogram(aes(fill = Diplotype), color = "black", bins = 100)+
ggeasy::easy_remove_legend_title()+
theme_bw(base_size = 14)+
xlab(bquote(V[Max]))+
ylab("Population (n)")
dist.grb <- dist.plot +
ggeasy::easy_remove_legend() +
scale_x_log10()+
xlab("Log Scale")+
ylab("Count (n)")
dist.plot2 <- ggplot(
data = population.sim(n=1000,`CV%_D` = 25, Start = 0.5,
By = 5, pop_freq = F) %>%
group_by(ID) %>%
filter(S == min(S)),
aes(x = Km))+
geom_histogram(aes(fill = Diplotype), color = "black", bins = 100)+
ggeasy::easy_remove_legend_title()+
theme_bw(base_size = 14)+
xlab(bquote(K[M](uM)))+
ylab("Population (n)")+
scale_x_log10()
dist.grb2 <- dist.plot +
ggeasy::easy_remove_legend() +
scale_x_log10()+
xlab("Log Scale")+
ylab("Count (n)")
## Vmax and Km Distribution Plots
ggpubr::ggarrange(dist.plot +
annotation_custom(ggplotGrob(dist.grb),
xmin = 1000, xmax = 2500,
ymin = 250, ymax = 1250),
dist.plot2, ncol = 1, common.legend = T, legend = "right")
Distribution of Michaelis-Menten Parameters for Virtual Population (n=9000) stratified by CYP2D6 genotype. Km represent the Michaelis-Menten constant, and Vmax represents maximal 4-OH Atomoxetine formation (pmol/min/mg protein) by CYP2D6 in human liver microsomes.
Virtual Population Vmax and Km Distribution Table
population.sim(n = 1000, End = 0,`CV%_D` = 25) %>%
group_by(Diplotype) %>%
summarise(Vmax_mean = round(mean(Vmax),2),
Vmax_sd = round(sd(Vmax),2),
Km_mean = round(mean(Km),2),
Km_sd = round(sd(Km),2),
`n` = length(Vmax),
Vmax = paste0(Vmax_mean," ± ",Vmax_sd),
Km = paste0(Km_mean," ± ",Km_sd)) %>%
ungroup() %>%
left_join(diplotype.data %>%
mutate(across(where(is.numeric), round, 2)) %>%
select(Diplotype, Vmax, Km) %>%
rename("Actual (Vmax)" = "Vmax", "Actual (Km)" = "Km"),
by = "Diplotype") %>%
select(Diplotype,`Actual (Vmax)`, Vmax, `Actual (Km)`,Km, n) %>%
rename("Simulated." = "Vmax",
"Simulated" = "Km",
"Actual." = "Actual (Vmax)",
"Actual" = "Actual (Km)") %>%
kbl(table.attr = "style='width:60%;'",
caption = "Simulated Population Parameter Estimates") %>%
kable_classic("striped", full_width = T) %>%
add_header_above(c(" ", "Vmax Estimate" = 2, "Km Estimate" = 2, " "))|
Vmax Estimate
|
Km Estimate
|
||||
|---|---|---|---|---|---|
| Diplotype | Actual. | Simulated. | Actual | Simulated | n |
| 1/1 | 354.50 | 352.17 ± 88.91 | 1.45 | 1.44 ± 0.36 | 1000 |
| 1/2 | 470.33 | 467.24 ± 117.97 | 1.39 | 1.38 ± 0.35 | 1000 |
| 1/2x2 | 1410.00 | 1400.73 ± 353.65 | 1.82 | 1.81 ± 0.46 | 1000 |
| 1/4 | 274.00 | 272.2 ± 68.72 | 0.92 | 0.91 ± 0.23 | 1000 |
| 2/2 | 252.00 | 250.34 ± 63.2 | 4.59 | 4.56 ± 1.15 | 1000 |
| 2/3 | 251.00 | 249.35 ± 62.95 | 1.32 | 1.31 ± 0.33 | 1000 |
| 2/4 | 95.90 | 95.27 ± 24.05 | 1.68 | 1.67 ± 0.42 | 1000 |
| 4/41 | 30.90 | 30.7 ± 7.75 | 5.35 | 5.31 ± 1.34 | 1000 |
| 4/5 | 12.32 | 12.24 ± 3.09 | 69.15 | 68.7 ± 17.34 | 1000 |
Incorporation of Experimental Error
Residual error was added to simulated data using a constant coefficient of variation structure, where the observed metabolic rate (Vobs)for the ithindividual at the jth substrate incubation concentration. The residual errors (εi) were normally distributed with a mean of 0 and a standard deviation (σ) set at the following values: 0, 0.05, 0.10, or 0.2; corresponding to a coefficient of variation (CV%) of 0%, 5%, 10%, and 20%.
# Error Design Datatable -------
error_design.data <- rbind(population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 0) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "0% CV"),
population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 5) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "5% CV"),
population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 10) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "10% CV"),
population.sim(n=3, Start = 0, By = 0.1, `CV%_V` = 20) %>%
filter(Diplotype == "1/1") %>% mutate(Condition = "20% CV")) %>%
filter(S %in% c(0.1, 0.5, 1, 2, 5, 10, 15, 25, 40, 55, 70, 85, 100)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Rich Design (9 pts)")
# Assigning Condition as Ordered Factor -----------
error_design.data$Condition <- factor(error_design.data$Condition,
levels = c("0% CV","5% CV","10% CV","20% CV"))
# Simulation Grid ---------
error.nls <- nls(formula = V ~ (Vmax*S)/(Km + S), data = error_design.data,
start = list(Vmax = 350, Km = 2))
error.grid <- tibble(S = seq(0,100,0.1),
V = (coef(error.nls)[[1]]*S)/(coef(error.nls)[[2]]+S))# Experimental Error Scenarios Plot -----------
ggplot(error_design.data %>% filter(S %in% c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)),
aes(x = S, y = V))+
geom_line(data = error.grid, size = 1, alpha = 0.75)+
geom_point(aes(y = V_adj, group = ID2, color = ID2, shape = ID2), size = 2.5)+
facet_wrap(~Condition, ncol = 2, scales = "free")+
theme_bw(base_size = 14)+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.7)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.9,0.11), legend.box = "horizontal",
legend.background = element_rect(color = "grey"))+
labs(title = "Experimental Error Scenarios")
Illustration of experimental error scenarios for 3 subjects with the same genotype. Black line represents the non-linear least-squares Michaelis-Menten fit for the genotype group.
In-Silico Experiments
Strategic Sampling Approach
A custom function, update_data(), was created to
generate experimental data in-silico according to the specified
experimental design and conditions. Two in-silico experimental designs
(rich design, and sparse design) for two experimental conditions (UM/EM
and IM/PM conditions) were incorporated into the function. Rich design
experiments were conducted using 9 overlapping atomoxetine (ATX)
incubation concentrations per subject, while sparse design experiments
used 3-4 staggered incubation concentrations per subject , with both
designs covering similar concentration ranges according to their
experimental condition (UM/EM range: 0.10 - 100 uM, IM/PM range: 1-2000
uM). Additionally, each experimental design can subjected to variable
degrees of experimental error using the CV% input. The
inputs for the function the update_data() function include
the following:
n_indiv - number of subjects per diplotype group (default = 10).
CV%_D - Average standard deviation for each diplotype group as a percentage of the true estimate for Vmax and Km (between subject variability for each diplotype group) (default = 25).
CV% - Average coefficient of variation across all points for simulated Michaelis-Menten kinetics for each individual (residual error).
start - Starting concentration from Michaelis-Menten simulation (default = 0).
end - Ending concentration for Michaelis-Menten kinetics (default = 100 uM).
by - Simulated concentration increment (default = 0.5).
Pop_freq - TRUE or FALSE argument indicating whether to use true population frequency (default is FALSE).
Seed - Random number generator seed (for reproducibility, default = 23457).
type - Experimental design to be simulated (options: rich sampling (9pt) = “rich”, sparse sampling (3pt) = “sparse3pt”, sparse sampling (4pt) = “sparse4pt”, rich sampling for poor metabolizers (16 pt) = “PM”, sparse sampling for poor metabolizer (3pt) = “PM_3pt”, sparse sampling for poor metabolizer (4pt) = “PM_4pt”.
update_data <-function(`CV%`, type, n_indiv = 10, `CV%.D` = 25, by = 0.1,
start = 0, end = 100, Seed = 23457, Pop.freq = FALSE){
# Sampling Information --------------
# Full Range Set
full_set <- c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)
PM_range <- c(1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, 2000)
# Strategic Sampling Sets ----
sparse3pt_set1 <- c(0.1, 2.5, 50)
sparse3pt_set2 <- c(1, 10, 100)
sparse4pt_set1 <- c(0.1, 2.5, 25, 50)
sparse4pt_set2 <- c(1, 10, 50, 100)
PM_3pt_set1 <- c(10,100,1000)
PM_3pt_set2 <- c(25,250,2000)
PM_4pt_set1 <- c(1,10,100,1000)
PM_4pt_set2 <- c(5,25,250,2000)
PM_range_set1 <- c(1,10,25,50,100,400,1000,2000)
PM_range_set2 <- c(5,15,30,60,90,200,800,1600)
# Rich Sampling Simulation ----------------
Rich <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
# Sparse Sampling Simulation (3 point) ------------
Sparse3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse3pt_set1, S %in% sparse3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
# Sparse Sampling Simulation (4 point) --------------
Sparse4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse4pt_set1, S %in% sparse4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
PM <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
PM_3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_3pt_set1, S %in% PM_3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
PM_4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_4pt_set1, S %in% PM_4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V)
# Output Selection ----------
switch(type,
"rich" = Rich,
"sparse3pt" = Sparse3pt,
"sparse4pt" = Sparse4pt,
"PM" = PM,
"PM_3pt" = PM_3pt,
"PM_4pt" = PM_4pt)
}Representative Experimental Design Plots
Strategic sampling was accomplished by evenly dividing individuals falling under the same CYP2D6 genotype into two groups (Group 1 and Group 2). Each group was designated a pre-defined set of incubation concentrations to be used for the in-silico experiment with the total range of incubation concentrations dependent on the type of CYP2D6 metabolizer (as described in the Strategic Sampling Approach section above). Note: strategic sampling was not used for rich sample designs, meaning all subject were subjected to the same set of incubations concentrations.
Extensive and Ultra-Rapid Metabolizers
Incubation Concentrations - Rich (9pt) - Set: 0.1, 0.5, 1, 2.5, 5, 10, 25, 50, and 100 uM
Incubation Concentrations - Sparse (3pt) - Set 1: 0.1, 2.5, and 50 uM - Set 2: 1, 10, and 100 uM
Incubation Concentrations - Sparse (4pt) - Set 1: 0.1, 2.5, 25, and 50 uM - Set 2: 1, 10, 50, and 100 uM
Intermediate and Poor Metabolizers
Incubation Concentrations - Rich (16pt) - Set: 1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, and 2000 uM
Incubation Concentrations - Sparse (3pt) - Set 1: 10,100, and 1000 uM - Set 2: 25, 250, and 2000 uM
Incubation Concentrations - Sparse (4pt) - Set 1: 1,10,100, and 1000 uM - Set 2: 5,25,250, and 2000 uM
# Generating Representative Data sets =========
## Representative data comes from 2 wild-type subjects.
rich_design.data <- update_data(`CV%` = 0, type = "rich") %>%
filter(Diplotype == "1/1", ID %in% c(1,2)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Rich Design (9 pts)")
sparse_3pt_design.data <- update_data(`CV%` = 0, type = "sparse3pt")%>%
filter(Diplotype == "1/1", ID %in% c(1,2)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (3 pts)")
sparse_4pt_design.data <- update_data(`CV%` = 0, type = "sparse4pt")%>%
filter(Diplotype == "1/1", ID %in% c(1,2)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (4 pts)")
PM_design.data <- update_data(`CV%` = 0, type = "PM") %>%
filter(Diplotype == "4/5", ID %in% c(81,82)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Rich Design (16 pts)")
PM3pt_design.data <- update_data(`CV%` = 0, type = "PM_3pt") %>%
filter(Diplotype == "4/5", ID %in% c(81,82)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (3 pts)")
PM4pt_design.data <- update_data(`CV%` = 0, type = "PM_4pt") %>%
filter(Diplotype == "4/5", ID %in% c(81,82)) %>%
mutate(ID2 = paste("Subject",ID),
`Study Design` = "Sparse Design (4 pts)")
# Solution Grid
design.nls <- nls(formula = V ~ (Vmax*S)/(Km + S), data = rich_design.data,
start = list(Vmax = 350, Km = 2))
grid <- tibble(S = seq(0,100,0.1),
V = (coef(design.nls)[[1]]*S)/(coef(design.nls)[[2]]+S))
design_PM.nls <- nls(formula = V ~ (Vmax*S)/(Km + S), data = PM_design.data,
start = list(Vmax = 20, Km = 50))
grid_pm <- tibble(S = seq(0,2000,0.1),
V = (coef(design_PM.nls)[[1]]*S)/(coef(design_PM.nls)[[2]]+S))# UM/EM STRATEGIC SAMPLING DESIGN FIGURE ================
## Representative plots contain data for 2 wild-type subjects.
# Rich Design--------
D1<- ggplot(rich_design.data, aes(x = S, y = V))+
geom_line(data = grid, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Conventional")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# Sparse Design (3pt)----------
D2<- ggplot(sparse_3pt_design.data, aes(x = S, y = V))+
geom_line(data = grid, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Strategic Scenario (2)")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# Sparse Design (4pt)------------
D3 <- ggplot(sparse_4pt_design.data, aes(x = S, y = V))+
geom_line(data = grid, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Strategic Scenario (1)")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# IM/PM STRATEGIC SAMPLING DESIGN FIGURE ================
## Representative plots contain data for 2 wild-type subjects.
# PM Sparse Design (4pt)---------
D4 <- ggplot(PM4pt_design.data, aes(x = S, y = V))+
geom_line(data = grid_pm, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# Extensize metabolizer model---------
D5 <- D3 + labs(title = "") #4pt
D6 <- D2 + labs(title = "") #3pt
D7 <- D1 + labs(title = "Extensive Metabolizer") #9pt
# PM Sparse Design (3pt)---------
D8 <- ggplot(PM3pt_design.data, aes(x = S, y = V))+
geom_line(data = grid_pm, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)
# PM Rich Design (16pt)---------
D9 <- ggplot(PM_design.data, aes(x = S, y = V))+
geom_line(data = grid_pm, size = 1.5, alpha = 0.75)+
geom_point(aes(color=ID2, shape = ID2), size = 3.5)+
theme_bw(base_size = 14)+
ggeasy::easy_remove_legend_title()+
theme(legend.position = c(0.7,0.2))+
facet_wrap(~`Study Design`, ncol = 1)+
labs(title = "Poor Metabolizers")+
xlab("[Substrate](uM)")+
ylab("Reaction Rate (V)\n")+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.6)UM/EM Design
## Extensive Metabolizers
ggpubr::ggarrange(D7,D5,D6,ncol = 2, nrow = 2, labels = "AUTO")
IM/PM Design
## Poor Metabolizers
ggpubr::ggarrange(D9,D4,D8, ncol = 2, nrow = 2, labels = "AUTO")
Population Modeling (UM/EM)
Generation of Experimental Data (Base Model - 0% CV)
# Population Datasets (CV = 0%)
rich.data <- update_data(`CV%` = 0, type = "rich") %>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt.data <- update_data(`CV%` = 0, type = "sparse3pt") %>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt.data <- update_data(`CV%` = 0, type = "sparse4pt") %>%
filter(!Diplotype %in% c("4/41", "4/5"))Individual-Fits of Experimental Data (Base Model - 0% CV)
A grouped data frame was created for each experiment where reaction
rate (V) is predicted by substrate concentration (S) according to
subject (ID). Individual fits (un-weighted) and parameter estimates for
the data sets were obtained using non-linear least squares function
nlsList(). Values obtained from the nlsList()
were used as initial estimates for the non-linear mixed effect
model.
# Grouped Data frame -------------
rich_data.grp <- groupedData(V~S|ID, rich.data)
sparse3pt.grp <- groupedData(V~S|ID, sparse3pt.data)
sparse4pt.grp <- groupedData(V~S|ID, sparse4pt.data)
# Fit Model ((Non-Linear Least Squares) ------------
rich_fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_data.grp)
sparse3pt_fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt.grp)
sparse4pt_fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt.grp)Non-Linear Mixed Effect Modeling (Base Model - 0% CV)
Non-linear mixed effect modeling was conducted using the
nlme() package version 3.1 (Bates et al., 2000). The
nlsList() fits were used as initial estimates for the NLME
model. Individual level random effects were applied to both Vmax and Km.
A diagonal variance-covariance structure was assigned to each model
using the pdDiag() function.
# Fit Model ((Non-Linear Mixed-Effect) ------------
rich.nlme <- nlme(rich_fit.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt.nlme <- nlme(sparse3pt_fit.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt.nlme <- nlme(sparse4pt_fit.nls, random = pdDiag(Vmax + Km ~ 1))
# Model Fixed Effect comparisons
fixef.table <- rbind(fixef(rich.nlme),
fixef(sparse3pt.nlme),
fixef(sparse4pt.nlme)) %>%
as_tibble() %>%
mutate(Model = c("Rich","Sparse (3pt)","Sparse (4pt)")) %>%
relocate(Model)
fixef.table %>%
mutate(across(where(is.numeric), round, 2)) %>%
kbl( table.attr = "style='width:60%;'",
caption = "Base Model: Mixed Effect Michaelis-Menten") %>%
kable_classic("striped")| Model | Vmax | Km |
|---|---|---|
| Rich | 459.63 | 1.95 |
| Sparse (3pt) | 459.63 | 1.95 |
| Sparse (4pt) | 459.63 | 1.95 |
Non-Linear Mixed Effect Modeling (Covariate Model - 0% CV)
CYP2D6 genotype was incorporated as a categorical covariate in the model by adding it as a fixed effect that impacts both Vmax and Km. The extracted population level fixed effects for Vmax and Km from the base model (described in the previous section) were used as initial estimates for the covariate model.
#Models with genotype as a Covariate -----
rich.fxf <- fixef.table[1,2:3]
sparse3pt.fxf <- fixef.table[2,2:3]
sparse4pt.fxf <- fixef.table[3,2:3]
rich_covar.nlme <- update(rich.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich.fxf[[1]], rep(0,6),
rich.fxf[[2]], rep(0,6)))
sparse3pt_covar.nlme <- update(sparse3pt.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt.fxf[[1]], rep(0,6),
sparse3pt.fxf[[2]], rep(0,6)))
sparse4pt_covar.nlme <- update(sparse4pt.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt.fxf[[1]], rep(0,6),
sparse4pt.fxf[[2]], rep(0,6)))Base Model and Covariate Model Evaluation (0% CV)
The impact of adding CYP2D6 genotype as a (covariate on both Vmax and
Km) and the whether adding covariates to the model improved fit compared
to baseline was assessed using analysis of variance
anova().
Overall Impact of adding covariates to base model
anova(rich.nlme, rich_covar.nlme)
anova(sparse3pt.nlme, sparse3pt_covar.nlme)
anova(sparse4pt.nlme, sparse4pt_covar.nlme)Impact of adding covariate to Vmax and Km
anova(rich_covar.nlme) %>%
kbl(caption = "Rich Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")
anova(sparse3pt_covar.nlme)%>%
kbl(caption = "Sparse 3pt Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")
anova(sparse4pt_covar.nlme)%>%
kbl(caption = "Sparse 4pt Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")Summary Tables of Fixed Effects Estimates (Covariate Model - 0% CV)
Table of Covariate Model Estimates
# Covariate Analysis Summary Table ------------
covaref.table <- rbind(fixef(rich_covar.nlme),
fixef(sparse3pt_covar.nlme),
fixef(sparse4pt_covar.nlme)) %>%
as_tibble() %>%
mutate(`Covariate Model` = c("Rich","Sparse (3pt)","Sparse (4pt)")) %>%
pivot_longer(`Vmax.(Intercept)`:`Km.Diplotype2/4`,
names_to = "Parameter",
values_to = "Value") %>%
pivot_wider(names_from = `Covariate Model`, values_from = Value) %>%
separate(Parameter, sep = "\\.", into = c("Parameter","Diplotype"), remove = T) %>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (1/1)"),
across(where(is.numeric), round, 2))Reference Group (*1/*1) Parameter Estimates
# Reference Estimates -----------
covaref.table %>%
filter(Diplotype == "Reference (1/1)")%>%
kbl(caption = "Reference Population Estimates",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")| Parameter | Diplotype | Rich | Sparse (3pt) | Sparse (4pt) |
|---|---|---|---|---|
| Vmax | Reference (1/1) | 367.01 | 367.01 | 367.01 |
| Km | Reference (1/1) | 1.50 | 1.50 | 1.50 |
Variant Group ETAs Table
# CYP2D6 Genotype Effects --------
covaref.table%>%
kbl(caption = "Michaelis Menten Mixed-Effect Model (w/Covariates)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped") %>%
pack_rows("Vmax Estimate (Population)", 1,1) %>%
pack_rows("Covariate (Vmax Eta)", 2, 7)%>%
pack_rows("Km Estimate (Population)", 8,8) %>%
pack_rows("Covariate (Km Eta)", 9, 14)| Parameter | Diplotype | Rich | Sparse (3pt) | Sparse (4pt) |
|---|---|---|---|---|
| Vmax Estimate (Population) | ||||
| Vmax | Reference (1/1) | 367.01 | 367.01 | 367.01 |
| Covariate (Vmax Eta) | ||||
| Vmax | Diplotype1/2 | 119.92 | 119.92 | 119.92 |
| Vmax | Diplotype1/2x2 | 1092.74 | 1092.74 | 1092.74 |
| Vmax | Diplotype1/4 | -83.34 | -83.34 | -83.34 |
| Vmax | Diplotype2/2 | -106.12 | -106.11 | -106.12 |
| Vmax | Diplotype2/3 | -107.15 | -107.15 | -107.15 |
| Vmax | Diplotype2/4 | -267.73 | -267.72 | -267.73 |
| Km Estimate (Population) | ||||
| Km | Reference (1/1) | 1.50 | 1.50 | 1.50 |
| Covariate (Km Eta) | ||||
| Km | Diplotype1/2 | -0.06 | -0.06 | -0.06 |
| Km | Diplotype1/2x2 | 0.39 | 0.39 | 0.39 |
| Km | Diplotype1/4 | -0.55 | -0.55 | -0.55 |
| Km | Diplotype2/2 | 3.25 | 3.25 | 3.25 |
| Km | Diplotype2/3 | -0.13 | -0.13 | -0.13 |
| Km | Diplotype2/4 | 0.23 | 0.23 | 0.23 |
Covariate Model 95% Confidence Interval Table (0% CV)
# Function to Extract Confidence Interval Info ----------
extract_CI <- function(int.data, name){
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(across(where(is.numeric), round, 2),
`CI (95%)` = paste0("(",lower,", ",upper,")")) %>%
select(est., `CI (95%)`)
colnames(temp) <- c(paste(name, " Estimate"), paste(name, " CI (95%)"))
temp
}
# Confidence Interval by Sampling Scenario -------
rich_covar.intervals <- intervals(rich_covar.nlme, which = "fixed")
sparse3pt_covar.intervals <- intervals(sparse3pt_covar.nlme, which = "fixed")
sparse4pt_covar.intervals <- intervals(sparse4pt_covar.nlme, which = "fixed")
# Summary Table of Confidence Intervals by Diplotype and Scenario -----------
covar_ci.table <- cbind(extract_CI(rich_covar.intervals, name = "Rich"),
extract_CI(sparse3pt_covar.intervals, name = "Sparse 3pt"),
extract_CI(sparse4pt_covar.intervals, name = "Sparse 4pt")) %>%
mutate(Parameter = rownames(rich_covar.intervals$fixed)) %>%
relocate(Parameter) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (1/1)"))
colnames(covar_ci.table)<- c("Parameter", "Diplotype",
"Estimate","CI (95%)",
"Estimate", "CI (95%)",
"Estimate", "CI (95%)")
covar_ci.table%>%
kbl(caption = "Michaelis Menten Mixed-Effect Model (w/Covariates)") %>%
kable_classic(full_width = T, "striped")%>%
pack_rows("Vmax Estimate (Population)", 1,1) %>%
pack_rows("Covariate (Vmax Eta)", 2, 7)%>%
pack_rows("Km Estimate (Population)", 8,8) %>%
pack_rows("Covariate (Km Eta)", 9, 14) %>%
add_header_above(c(" ", " ",
"Rich (9pt)" = 2,
"Sparse (3pt)" = 2,
"Sparse (4pt)" = 2))|
Rich (9pt)
|
Sparse (3pt)
|
Sparse (4pt)
|
|||||
|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Estimate | CI (95%) | Estimate | CI (95%) | Estimate | CI (95%) |
| Vmax Estimate (Population) | |||||||
| Vmax | Reference (1/1) | 367.01 | (319.69, 414.33) | 367.01 | (319.34, 414.67) | 367.01 | (319.51, 414.51) |
| Covariate (Vmax Eta) | |||||||
| Vmax | Diplotype1/2 | 119.92 | (53.01, 186.84) | 119.92 | (52.51, 187.33) | 119.92 | (52.74, 187.1) |
| Vmax | Diplotype1/2x2 | 1092.74 | (1025.83, 1159.66) | 1092.74 | (1025.33, 1160.16) | 1092.74 | (1025.56, 1159.93) |
| Vmax | Diplotype1/4 | -83.34 | (-150.26, -16.42) | -83.34 | (-150.75, -15.93) | -83.34 | (-150.52, -16.16) |
| Vmax | Diplotype2/2 | -106.12 | (-173.03, -39.2) | -106.11 | (-173.52, -38.7) | -106.12 | (-173.3, -38.94) |
| Vmax | Diplotype2/3 | -107.15 | (-174.07, -40.24) | -107.15 | (-174.56, -39.74) | -107.15 | (-174.33, -39.97) |
| Vmax | Diplotype2/4 | -267.73 | (-334.64, -200.81) | -267.72 | (-335.13, -200.31) | -267.73 | (-334.91, -200.54) |
| Km Estimate (Population) | |||||||
| Km | Reference (1/1) | 1.50 | (1.33, 1.67) | 1.50 | (1.33, 1.67) | 1.50 | (1.33, 1.67) |
| Covariate (Km Eta) | |||||||
| Km | Diplotype1/2 | -0.06 | (-0.31, 0.18) | -0.06 | (-0.31, 0.18) | -0.06 | (-0.31, 0.18) |
| Km | Diplotype1/2x2 | 0.39 | (0.14, 0.63) | 0.39 | (0.14, 0.63) | 0.39 | (0.14, 0.63) |
| Km | Diplotype1/4 | -0.55 | (-0.79, -0.3) | -0.55 | (-0.79, -0.3) | -0.55 | (-0.79, -0.3) |
| Km | Diplotype2/2 | 3.25 | (3.01, 3.49) | 3.25 | (3.01, 3.5) | 3.25 | (3.01, 3.5) |
| Km | Diplotype2/3 | -0.13 | (-0.38, 0.11) | -0.13 | (-0.38, 0.11) | -0.13 | (-0.38, 0.11) |
| Km | Diplotype2/4 | 0.23 | (-0.01, 0.48) | 0.23 | (-0.01, 0.48) | 0.23 | (-0.01, 0.48) |
Impact of Residual Error
Data sets which include increasing degrees of residual error (5, 10,
and 20% CV) were generated using the update_data() function
described earlier (see In-Silico Experiments section for full
description). Additonally, the same methodology for model analysis used
in the the previous section (Population Michaelis-Menten Modeling
(UM/EM)) was applied here to evaluate the impact of increasing
experimental error on parameter estimates.
Generation of Experimental Data (Base Model - 5-20% CV)
# Population Data sets (CV = 5%) ----------------
rich_CV5.data <- update_data(`CV%` = 5, type = "rich")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt_CV5.data <- update_data(`CV%` = 5, type = "sparse3pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt_CV5.data <- update_data(`CV%` = 5, type = "sparse4pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
# Population Data sets (CV = 10%) ----------------
rich_CV10.data <- update_data(`CV%` = 10, type = "rich")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt_CV10.data <- update_data(`CV%` = 10, type = "sparse3pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt_CV10.data <- update_data(`CV%` = 10, type = "sparse4pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
# Population Data sets (CV = 20%) ----------------
rich_CV20.data <- update_data(`CV%` = 20, type = "rich")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse3pt_CV20.data <- update_data(`CV%` = 20, type = "sparse3pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))
sparse4pt_CV20.data <- update_data(`CV%` = 20, type = "sparse4pt")%>%
filter(!Diplotype %in% c("4/41", "4/5"))Inidividual-Fits of Experimental Data (Base Model - 5-20% CV)
# Grouped Data frame =============================
rich_CV5_data.grp <- groupedData(V~S|ID, rich_CV5.data)
sparse3pt_CV5.grp <- groupedData(V~S|ID, sparse3pt_CV5.data)
sparse4pt_CV5.grp <- groupedData(V~S|ID, sparse4pt_CV5.data)
rich_CV10_data.grp <- groupedData(V~S|ID, rich_CV10.data)
sparse3pt_CV10.grp <- groupedData(V~S|ID, sparse3pt_CV10.data)
sparse4pt_CV10.grp <- groupedData(V~S|ID, sparse4pt_CV10.data)
rich_CV20_data.grp <- groupedData(V~S|ID, rich_CV20.data)
sparse3pt_CV20.grp <- groupedData(V~S|ID, sparse3pt_CV20.data)
sparse4pt_CV20.grp <- groupedData(V~S|ID, sparse4pt_CV20.data)
# Fit Model ((Non-Linear Least Squares) =====================
rich_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_CV5_data.grp)
sparse3pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt_CV5.grp)
sparse4pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt_CV5.grp)
rich_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_CV10_data.grp)
sparse3pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt_CV10.grp)
sparse4pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt_CV10.grp)
rich_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = rich_CV20_data.grp)
sparse3pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse3pt_CV20.grp)
sparse4pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = sparse4pt_CV20.grp)Non-Linear Mixed Effect Modeling (Base Model - 5-20% CV)
# Fit Model ((Non-Linear Mixed-Effect) ===================
rich_CV5.nlme <- nlme(rich_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt_CV5.nlme <- nlme(sparse3pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt_CV5.nlme <- nlme(sparse4pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
rich_CV10.nlme <- nlme(rich_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt_CV10.nlme <- nlme(sparse3pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt_CV10.nlme <- nlme(sparse4pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
rich_CV20.nlme <- nlme(rich_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
sparse3pt_CV20.nlme <- nlme(sparse3pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
sparse4pt_CV20.nlme <- nlme(sparse4pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))Impact of Residual Error on Fixed Effects
Summary Table of Fixed Effects (w/Residual Error)
# Model Fixed Effect comparisons----------
fixef_CV.table <- rbind(
fixef(rich.nlme),
fixef(sparse3pt.nlme),
fixef(sparse4pt.nlme),
fixef(rich_CV5.nlme),
fixef(sparse3pt_CV5.nlme),
fixef(sparse4pt_CV5.nlme),
fixef(rich_CV10.nlme),
fixef(sparse3pt_CV10.nlme),
fixef(sparse4pt_CV10.nlme),
fixef(rich_CV20.nlme),
fixef(sparse3pt_CV20.nlme),
fixef(sparse4pt_CV20.nlme)) %>%
as_tibble() %>%
mutate(Model = c("Rich (0% CV)","Sparse (0% CV, 3pt)","Sparse (0% CV, 4pt)",
"Rich (5% CV)","Sparse (5% CV, 3pt)","Sparse (5% CV, 4pt)",
"Rich (10% CV)","Sparse (10% CV, 3pt)","Sparse (10% CV, 4pt)",
"Rich (20% CV)","Sparse (20% CV, 3pt)","Sparse (20% CV, 4pt)"),
across(where(is.numeric),round,3)) %>%
relocate(Model)
# Fixed Effect Table across conditions -------------------
fixef_CV.table %>%
kbl(caption = "Population Estimates (w/ Residual Error)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")%>%
pack_rows(" ", 1,3) %>%
pack_rows(" ", 4, 6)%>%
pack_rows(" ", 7, 9)%>%
pack_rows(" ", 10, 12)| Model | Vmax | Km |
|---|---|---|
| Rich (0% CV) | 459.628 | 1.948 |
| Sparse (0% CV, 3pt) | 459.628 | 1.948 |
| Sparse (0% CV, 4pt) | 459.628 | 1.948 |
| Rich (5% CV) | 458.428 | 1.689 |
| Sparse (5% CV, 3pt) | 460.684 | 1.726 |
| Sparse (5% CV, 4pt) | 459.051 | 1.692 |
| Rich (10% CV) | 461.150 | 1.723 |
| Sparse (10% CV, 3pt) | 465.169 | 1.888 |
| Sparse (10% CV, 4pt) | 461.482 | 1.804 |
| Rich (20% CV) | 464.268 | 1.765 |
| Sparse (20% CV, 3pt) | 468.840 | 1.924 |
| Sparse (20% CV, 4pt) | 464.632 | 1.924 |
Summary Table of Fixed Effects (w/Residual Error & 95% CI)
# Function to extract confidence intervals from model ---------
extract_CI2 <- function(int.data, name){
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(across(where(is.numeric), round, 2),
`CI (95%)` = paste0("(",lower,", ",upper,")")) %>%
select(est., `CI (95%)`)
colnames(temp) <- c("Estimate","CI (95%)")
temp$Model <- name
temp$Parameter <- c("Vmax","Km")
temp
}
# Confidence Interval Data ---------
rich_CV5.intervals <- intervals(rich_CV5.nlme, which = "fixed")
sparse3pt_CV5.intervals <- intervals(sparse3pt_CV5.nlme, which = "fixed")
sparse4pt_CV5.intervals <- intervals(sparse4pt_CV5.nlme, which = "fixed")
rich_CV10.intervals <- intervals(rich_CV10.nlme, which = "fixed")
sparse3pt_CV10.intervals <- intervals(sparse3pt_CV10.nlme, which = "fixed")
sparse4pt_CV10.intervals <- intervals(sparse4pt_CV10.nlme, which = "fixed")
rich_CV20.intervals <- intervals(rich_CV20.nlme, which = "fixed")
sparse3pt_CV20.intervals <- intervals(sparse3pt_CV20.nlme, which = "fixed")
sparse4pt_CV20.intervals <- intervals(sparse4pt_CV20.nlme, which = "fixed")
# Confidence Interval Tables ----------
CV5_ci.table <- rbind(extract_CI2(rich_CV5.intervals, name = "Rich (5% CV)"),
extract_CI2(sparse3pt_CV5.intervals, name = "Sparse 3pt (5% CV)"),
extract_CI2(sparse4pt_CV5.intervals, name = "Sparse 4pt (5% CV)")) %>%
relocate(Model, Parameter) %>%
arrange(desc(Parameter))
CV10_ci.table <- rbind(extract_CI2(rich_CV10.intervals, name = "Rich (10% CV)"),
extract_CI2(sparse3pt_CV10.intervals, name = "Sparse 3pt (10% CV)"),
extract_CI2(sparse4pt_CV10.intervals, name = "Sparse 4pt (10% CV)")) %>%
relocate(Model, Parameter) %>%
arrange(desc(Parameter))
CV20_ci.table <- rbind(extract_CI2(rich_CV20.intervals, name = "Rich (20% CV)"),
extract_CI2(sparse3pt_CV20.intervals, name = "Sparse 3pt (20% CV)"),
extract_CI2(sparse4pt_CV20.intervals, name = "Sparse 4pt (20% CV)")) %>%
relocate(Model, Parameter) %>%
arrange(desc(Parameter))
# Summary Table of Estimates with CIs ----------
rbind(CV5_ci.table, CV10_ci.table, CV20_ci.table) %>%
arrange(desc(Parameter))%>%
kbl(caption = "Population Estimates (w/ Residual Error & 95% CI)",
table.attr = "style='width:60%;'") %>%
kable_classic(full_width = T, "striped")%>%
pack_rows("Vmax Estimates", 1,9) %>%
pack_rows("", 1, 3) %>%
pack_rows("", 4, 6) %>%
pack_rows("", 7, 9) %>%
pack_rows("Km Estimates", 10, 18) %>%
pack_rows(" ", 10, 12) %>%
pack_rows(" ", 13, 15) %>%
pack_rows(" ", 16, 18)| Model | Parameter | Estimate | CI (95%) |
|---|---|---|---|
| Vmax Estimates | |||
| Rich (5% CV) | Vmax | 458.43 | (356.53, 560.33) |
| Sparse 3pt (5% CV) | Vmax | 460.68 | (358.37, 563) |
| Sparse 4pt (5% CV) | Vmax | 459.05 | (357.03, 561.07) |
| Rich (10% CV) | Vmax | 461.15 | (358.84, 563.46) |
| Sparse 3pt (10% CV) | Vmax | 465.17 | (362.37, 567.97) |
| Sparse 4pt (10% CV) | Vmax | 461.48 | (358.88, 564.08) |
| Rich (20% CV) | Vmax | 464.27 | (360.44, 568.1) |
| Sparse 3pt (20% CV) | Vmax | 468.84 | (363.07, 574.61) |
| Sparse 4pt (20% CV) | Vmax | 464.63 | (360.16, 569.1) |
| Km Estimates | |||
| Rich (5% CV) | Km | 1.69 | (1.52, 1.86) |
| Sparse 3pt (5% CV) | Km | 1.73 | (1.59, 1.86) |
| Sparse 4pt (5% CV) | Km | 1.69 | (1.55, 1.83) |
| Rich (10% CV) | Km | 1.72 | (1.59, 1.85) |
| Sparse 3pt (10% CV) | Km | 1.89 | (1.74, 2.03) |
| Sparse 4pt (10% CV) | Km | 1.80 | (1.62, 1.99) |
| Rich (20% CV) | Km | 1.77 | (1.61, 1.93) |
| Sparse 3pt (20% CV) | Km | 1.92 | (1.66, 2.19) |
| Sparse 4pt (20% CV) | Km | 1.92 | (1.66, 2.19) |
Non-Linear Mixed Effect Modeling (Covariate Model - 5-20% CV)
# Models with Diplotype as a Covariate -----
rich_CV5.fxf <- fixef_CV.table[1,2:3]
sparse3pt_CV5.fxf <- fixef_CV.table[2,2:3]
sparse4pt_CV5.fxf <- fixef_CV.table[3,2:3]
rich_CV10.fxf <- fixef_CV.table[4,2:3]
sparse3pt_CV10.fxf <- fixef_CV.table[5,2:3]
sparse4pt_CV10.fxf <- fixef_CV.table[6,2:3]
rich_CV20.fxf <- fixef_CV.table[7,2:3]
sparse3pt_CV20.fxf <- fixef_CV.table[8,2:3]
sparse4pt_CV20.fxf <- fixef_CV.table[9,2:3]
# NLME Covariate Model w/5% CV ==================
rich_CV5_covar.nlme <- update(rich_CV5.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich_CV5.fxf[[1]], rep(0,6),
rich_CV5.fxf[[2]], rep(0,6)))
sparse3pt_CV5_covar.nlme <- update(sparse3pt_CV5.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt_CV5.fxf[[1]], rep(0,6),
sparse3pt_CV5.fxf[[2]], rep(0,6)))
sparse4pt_CV5_covar.nlme <- update(sparse4pt_CV5.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt_CV5.fxf[[1]], rep(0,6),
sparse4pt_CV5.fxf[[2]], rep(0,6)))
# NLME Covariate Model w/10% CV ==================
rich_CV10_covar.nlme <- update(rich_CV10.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich_CV10.fxf[[1]], rep(0,6),
rich_CV10.fxf[[2]], rep(0,6)))
sparse3pt_CV10_covar.nlme <- update(sparse3pt_CV10.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt_CV10.fxf[[1]], rep(0,6),
sparse3pt_CV10.fxf[[2]], rep(0,6)))
sparse4pt_CV10_covar.nlme <- update(sparse4pt_CV10.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt_CV10.fxf[[1]], rep(0,6),
sparse4pt_CV10.fxf[[2]], rep(0,6)))
# NLME Covariate Model w/20% CV ==================
rich_CV20_covar.nlme <- update(rich_CV20.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(rich_CV20.fxf[[1]], rep(0,6),
rich_CV20.fxf[[2]], rep(0,6)))
sparse3pt_CV20_covar.nlme <- update(sparse3pt_CV20.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse3pt_CV20.fxf[[1]], rep(0,6),
sparse3pt_CV20.fxf[[2]], rep(0,6)))
sparse4pt_CV20_covar.nlme <- update(sparse4pt_CV20.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(sparse4pt_CV20.fxf[[1]], rep(0,6),
sparse4pt_CV20.fxf[[2]], rep(0,6)))Covariate Models Evaluation (5-20% CV)
Covariate Model (5% CV)
# Summary of Covariate Models (w/ variability) -------
anova(rich_CV5_covar.nlme) # RICH SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 547 2321.6898 <.0001
## Vmax.Diplotype 6 547 370.0573 <.0001
## Km.(Intercept) 1 547 2549.1590 <.0001
## Km.Diplotype 6 547 30.1046 <.0001anova(sparse3pt_CV5_covar.nlme) # SPARSE 3PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 127 2590.3045 <.0001
## Vmax.Diplotype 6 127 366.3843 <.0001
## Km.(Intercept) 1 127 2536.6700 <.0001
## Km.Diplotype 6 127 19.4802 <.0001anova(sparse4pt_CV5_covar.nlme) # SPARSE 4PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 197 2515.2320 <.0001
## Vmax.Diplotype 6 197 357.4567 <.0001
## Km.(Intercept) 1 197 2596.3871 <.0001
## Km.Diplotype 6 197 19.9836 <.0001Covariate Model (10% CV)
anova(rich_CV10_covar.nlme) # RICH SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 547 2171.6492 <.0001
## Vmax.Diplotype 6 547 310.4159 <.0001
## Km.(Intercept) 1 547 1763.7233 <.0001
## Km.Diplotype 6 547 11.3821 <.0001anova(sparse3pt_CV10_covar.nlme) # SPARSE 3PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 127 2382.5424 <.0001
## Vmax.Diplotype 6 127 340.1663 <.0001
## Km.(Intercept) 1 127 756.2901 <.0001
## Km.Diplotype 6 127 6.9697 <.0001anova(sparse4pt_CV10_covar.nlme) # SPARSE 4PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 197 2366.5954 <.0001
## Vmax.Diplotype 6 197 339.9690 <.0001
## Km.(Intercept) 1 197 741.8162 <.0001
## Km.Diplotype 6 197 6.7521 <.0001Covariate Model (20% CV)
anova(rich_CV20_covar.nlme) # RICH SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 547 1768.6179 <.0001
## Vmax.Diplotype 6 547 258.1663 <.0001
## Km.(Intercept) 1 547 457.3300 <.0001
## Km.Diplotype 6 547 2.9774 0.0072anova(sparse3pt_CV20_covar.nlme) # SPARSE 3PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 127 1490.6401 <.0001
## Vmax.Diplotype 6 127 217.4530 <.0001
## Km.(Intercept) 1 127 193.1641 <.0001
## Km.Diplotype 6 127 2.7206 0.0161anova(sparse4pt_CV20_covar.nlme) # SPARSE 4PT SAMPLING## numDF denDF F-value p-value
## Vmax.(Intercept) 1 197 1716.8079 <.0001
## Vmax.Diplotype 6 197 252.6265 <.0001
## Km.(Intercept) 1 197 189.6812 <.0001
## Km.Diplotype 6 197 2.5226 0.0225Weighting Model (UM/EM)
Evaluating Heteroscedacity of Residual Error (UM/EM)
Heterogeneity of variance was examined qualitatively by plotting
standardized residules (Pearson). Heteroscedacity was corrected by
modeleling residual variance as a power function of mean fitted value
using the varPower() function.
Rich (0% CV)
plot(rich_covar.nlme, main = "Rich (9pt, 0% CV)")
plot(rich_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 0% CV, weighted)")
Sparse (3pt, 0% CV)
plot(sparse3pt_covar.nlme, main = "Sparse (3pt, 0% CV)")
plot(sparse3pt_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 0% CV, weighted)")
Sparse (4pt, 0% CV)
plot(sparse4pt_covar.nlme, main = "Sparse (4pt, 0% CV)")
plot(sparse4pt_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 0% CV, weighted)")
Rich (5% CV)
plot(rich_CV5_covar.nlme, main = "Rich (9pt, 5% CV)")
plot(rich_CV5_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 5% CV, weighted)")
Sparse (3pt, 5% CV)
plot(sparse3pt_CV5_covar.nlme, main = "Sparse (3pt, 5% CV)")
plot(sparse3pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 5% CV, weighted)")
Sparse (4pt, 5% CV)
plot(sparse4pt_CV5_covar.nlme, main = "Sparse (4pt, 5% CV)")
plot(sparse4pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 5% CV, weighted)")
Rich (10% CV)
plot(rich_CV10_covar.nlme, main = "Rich (9pt, 10% CV)")
plot(rich_CV10_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 10% CV, weighted)")
Sparse (3pt, 10% CV)
plot(sparse3pt_CV10_covar.nlme, main = "Sparse (3pt, 10% CV)")
plot(sparse3pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 10% CV, weighted)")
Sparse (4pt, 10% CV)
plot(sparse4pt_CV10_covar.nlme, main = "Sparse (4pt, 10% CV)")
plot(sparse4pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 10% CV, weighted)")
Rich (20% CV)
plot(rich_CV20_covar.nlme, main = "Rich (9pt, 20% CV)")
plot(rich_CV20_covar.nlme %>% update(weights = varPower()),
main = "Rich (9pt, 20% CV, weighted)")
Sparse (3pt, 20% CV)
plot(sparse3pt_CV20_covar.nlme, main = "Sparse (3pt, 20% CV)")
plot(sparse3pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "Sparse (3pt, 20% CV, weighted)")
Sparse (4pt, 20% CV)
plot(sparse4pt_CV20_covar.nlme, main = "Sparse (4pt, 20% CV)")
plot(sparse4pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "Sparse (4pt, 20% CV, weighted)")
Adjusting Residual Error Models
# Adjusting Residual Error Model ----
rich_CV5_covar.nlme <- rich_CV5_covar.nlme %>% update(weights = varPower())
sparse3pt_CV5_covar.nlme <- sparse3pt_CV5_covar.nlme %>% update(weights = varPower())
sparse4pt_CV5_covar.nlme <- sparse4pt_CV5_covar.nlme %>% update(weights = varPower())
rich_CV10_covar.nlme <- rich_CV10_covar.nlme %>% update(weights = varPower())
sparse3pt_CV10_covar.nlme <- sparse3pt_CV10_covar.nlme %>% update(weights = varPower())
sparse4pt_CV10_covar.nlme <- sparse4pt_CV10_covar.nlme %>% update(weights = varPower())
rich_CV20_covar.nlme <- rich_CV20_covar.nlme %>% update(weights = varPower())
sparse3pt_CV20_covar.nlme <- sparse3pt_CV20_covar.nlme %>% update(weights = varPower())
sparse4pt_CV20_covar.nlme <- sparse4pt_CV20_covar.nlme %>% update(weights = varPower())Covariate Model Summary Data
# Creating covariate fixed effect table ------
covaref_CV.table <- rbind(fixef(rich_CV5_covar.nlme),
fixef(sparse3pt_CV5_covar.nlme),
fixef(sparse4pt_CV5_covar.nlme),
fixef(rich_CV10_covar.nlme),
fixef(sparse3pt_CV10_covar.nlme),
fixef(sparse4pt_CV10_covar.nlme),
fixef(rich_CV20_covar.nlme),
fixef(sparse3pt_CV20_covar.nlme),
fixef(sparse4pt_CV20_covar.nlme)) %>%
as_tibble() %>%
mutate(`Covariate Model` = c("Rich (5% CV)","Sparse (5% CV, 3pt)","Sparse (5% CV, 4pt)",
"Rich (10% CV)","Sparse (10% CV, 3pt)","Sparse (10% CV, 4pt)",
"Rich (20% CV)","Sparse (20% CV, 3pt)","Sparse (20% CV, 4pt)")) %>%
pivot_longer(`Vmax.(Intercept)`:`Km.Diplotype2/4`,
names_to = "Parameter",
values_to = "Value") %>%
pivot_wider(names_from = `Covariate Model`, values_from = Value) %>%
separate(Parameter, sep = "\\.", into = c("Parameter","Diplotype"), remove = T) %>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (1/1)"),
across(where(is.numeric), round, 2))
# Complete Covariate Effects Table --------------------
covar_Vmax.table <- cbind(covaref.table, covaref_CV.table %>%
select(-Parameter, -Diplotype)) %>% filter(Parameter == "Vmax")
covar_Km.table <- cbind(covaref.table, covaref_CV.table %>%
select(-Parameter, -Diplotype)) %>% filter(Parameter == "Km")
# Vmax Estimates -------------------------------
final_vmax.est <- covar_Vmax.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (1/1)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (1/1)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (1/1)" = "Diplotype1/1")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Km), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Vmax)/Vmax)*100, digits = 2),
Vmax = round(Vmax, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Vmax, Predicted, `Error (%)`) %>%
rename("Predicted (Vmax)" = "Predicted") %>%
arrange(Condition)
# Km Estimates -----------------------------------
final_km.est <- covar_Km.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (1/1)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (1/1)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (1/1)" = "Diplotype1/1")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Vmax), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Km)/Km)*100, digits = 2),
Km = round(Km, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Km, Predicted, `Error (%)`) %>%
rename("Predicted (Km)" = "Predicted") %>%
arrange(Condition)Population Modeling (IM/PM)
Non-Linear Mixed Effect Modeling (Covariate Model - 0% CV)
# Full Population Simulation rich at 0% CV ------------
# Rich Sampling (16 pt, PM) ---------------
PM.data <- update_data(n_indiv = 10, `CV%` = 0, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM.grp <- groupedData(V~S|ID, PM.data)
PM.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM.grp)
PM.nlme <- nlme(PM.nls, random = pdDiag(Vmax + Km ~ 1))
PM_covar.nlme <- update(PM.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM.nlme)[[1]], rep(0,1),
fixed.effects(PM.nlme)[[2]], rep(0,1)))
# Sparse Sampling (3pt, PM) -------------
PM_3pt.data <- update_data(n_indiv = 10, `CV%` = 0, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt.grp <- groupedData(V~S|ID, PM_3pt.data)
PM_3pt.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt.grp)
PM_3pt.nlme <- nlme(PM_3pt.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_covar.nlme <- update(PM_3pt.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt.nlme)[[2]], rep(0,1)))
# Sparse Sampling (4pt, PM) ---------------
PM_4pt.data <- update_data(n_indiv = 10, `CV%` = 0, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt.grp <- groupedData(V~S|ID, PM_4pt.data)
PM_4pt.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt.grp)
PM_4pt.nlme <- nlme(PM_4pt.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_covar.nlme <- update(PM_4pt.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt.nlme)[[2]], rep(0,1)))Non-Linear Mixed Effect Modeling (Covariate Model - 5-20% CV)
# 5% Residual Error ---------------------------
PM_CV5.data <- update_data(n_indiv = 10, `CV%` = 5, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_CV5_data.grp <- groupedData(V~S|ID, PM_CV5.data)
PM_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_CV5_data.grp)
PM_CV5.nlme <- nlme(PM_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
PM_CV5_covar.nlme <- update(PM_CV5.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_CV5.nlme)[[1]], rep(0,1),
fixed.effects(PM_CV5.nlme)[[2]], rep(0,1)))
PM_3pt_CV5.data <- update_data(n_indiv = 10, `CV%` = 5, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt_CV5_data.grp <- groupedData(V~S|ID, PM_3pt_CV5.data)
PM_3pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt_CV5_data.grp)
PM_3pt_CV5.nlme <- nlme(PM_3pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_CV5_covar.nlme <- update(PM_3pt_CV5.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt_CV5.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt_CV5.nlme)[[2]], rep(0,1)))
PM_4pt_CV5.data <- update_data(n_indiv = 10, `CV%` = 5, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt_CV5_data.grp <- groupedData(V~S|ID, PM_4pt_CV5.data)
PM_4pt_CV5.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt_CV5_data.grp)
PM_4pt_CV5.nlme <- nlme(PM_4pt_CV5.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_CV5_covar.nlme <- update(PM_4pt_CV5.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt_CV5.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt_CV5.nlme)[[2]], rep(0,1)))
# 10% Residual Error --------------------
PM_CV10.data <- update_data(n_indiv = 10, `CV%` = 10, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_CV10_data.grp <- groupedData(V~S|ID, PM_CV10.data)
PM_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_CV10_data.grp)
PM_CV10.nlme <- nlme(PM_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
PM_CV10_covar.nlme <- update(PM_CV10.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_CV10.nlme)[[1]], rep(0,1),
fixed.effects(PM_CV10.nlme)[[2]], rep(0,1)))
PM_3pt_CV10.data <- update_data(n_indiv = 10, `CV%` = 10, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt_CV10_data.grp <- groupedData(V~S|ID, PM_3pt_CV10.data)
PM_3pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt_CV10_data.grp)
PM_3pt_CV10.nlme <- nlme(PM_3pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_CV10_covar.nlme <- update(PM_3pt_CV10.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt_CV10.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt_CV10.nlme)[[2]], rep(0,1)))
PM_4pt_CV10.data <- update_data(n_indiv = 10, `CV%` = 10, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt_CV10_data.grp <- groupedData(V~S|ID, PM_4pt_CV10.data)
PM_4pt_CV10.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt_CV10_data.grp)
PM_4pt_CV10.nlme <- nlme(PM_4pt_CV10.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_CV10_covar.nlme <- update(PM_4pt_CV10.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt_CV10.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt_CV10.nlme)[[2]], rep(0,1)))
# 20% Residual Error ----------
PM_CV20.data <- update_data(n_indiv = 10, `CV%` = 20, type = "PM") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_CV20_data.grp <- groupedData(V~S|ID, PM_CV20.data)
PM_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_CV20_data.grp)
PM_CV20.nlme <- nlme(PM_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
PM_CV20_covar.nlme <- update(PM_CV20.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_CV20.nlme)[[1]], rep(0,1),
fixed.effects(PM_CV20.nlme)[[2]], rep(0,1)))
PM_3pt_CV20.data <- update_data(n_indiv = 10, `CV%` = 20, type = "PM_3pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_3pt_CV20_data.grp <- groupedData(V~S|ID, PM_3pt_CV20.data)
PM_3pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_3pt_CV20_data.grp)
PM_3pt_CV20.nlme <- nlme(PM_3pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
PM_3pt_CV20_covar.nlme <- update(PM_3pt_CV20.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_3pt_CV20.nlme)[[1]], rep(0,1),
fixed.effects(PM_3pt_CV20.nlme)[[2]], rep(0,1)))
PM_4pt_CV20.data <- update_data(n_indiv = 10, `CV%` = 20, type = "PM_4pt") %>%
filter(Diplotype %in% c("4/41", "4/5"))
PM_4pt_CV20_data.grp <- groupedData(V~S|ID, PM_4pt_CV20.data)
PM_4pt_CV20.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = PM_4pt_CV20_data.grp)
PM_4pt_CV20.nlme <- nlme(PM_4pt_CV20.nls, random = pdDiag(Vmax + Km ~ 1))
PM_4pt_CV20_covar.nlme <- update(PM_4pt_CV20.nlme,
fixed = Vmax + Km ~ Diplotype,
start = c(fixed.effects(PM_4pt_CV20.nlme)[[1]], rep(0,1),
fixed.effects(PM_4pt_CV20.nlme)[[2]], rep(0,1)))Weighting Model (IM/PM)
Evaluating Heteroscedacity of Residual Error (IM/PM)
Rich (0% CV)
plot(PM_covar.nlme, main = "PM Rich (16pt, 0% CV)")
plot(PM_covar.nlme %>% update(weights = varPower()),
main = "PM Rich (16pt, 0% CV, weighted)")
Sparse (3pt, 0% CV)
plot(PM_3pt_covar.nlme, main = "PM Sparse (3pt, 0% CV)")
plot(PM_3pt_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 0% CV, weighted)")
Sparse (4pt, 0% CV)
plot(PM_4pt_covar.nlme, main = "PM Sparse (4pt, 0% CV)")
plot(PM_4pt_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 0% CV, weighted)")
Rich (5% CV)
plot(PM_CV5_covar.nlme, main = "PM Rich (16pt, 5% CV)")
plot(PM_CV5_covar.nlme %>% update(weights = varPower()),
main = "Rich (16pt, 5% CV, weighted)")
Sparse (3pt, 5% CV)
plot(PM_3pt_CV5_covar.nlme, main = "PM Sparse (3pt, 5% CV)")
plot(PM_3pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 5% CV, weighted)")
Sparse (4pt, 5% CV)
plot(PM_4pt_CV5_covar.nlme, main = "PM Sparse (4pt, 5% CV)")
plot(PM_4pt_CV5_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 5% CV, weighted)")
Rich (10% CV)
plot(PM_CV10_covar.nlme, main = "PM Rich (16pt, 10% CV)")
plot(PM_CV10_covar.nlme %>% update(weights = varPower()),
main = "PM Rich (16pt, 10% CV, weighted)")
Sparse (3pt, 10% CV)
plot(PM_3pt_CV10_covar.nlme, main = "PM Sparse (3pt, 10% CV)")
plot(PM_3pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 10% CV, weighted)")
Sparse (4pt, 10% CV)
plot(PM_4pt_CV10_covar.nlme, main = "PM Sparse (4pt, 10% CV)")
plot(PM_4pt_CV10_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 10% CV, weighted)")
Rich (20% CV)
plot(PM_CV20_covar.nlme, main = "PM Rich (16pt, 20% CV)")
plot(PM_CV20_covar.nlme %>% update(weights = varPower()),
main = "PM Rich (16pt, 20% CV, weighted)")
Sparse (3pt, 20% CV)
plot(PM_3pt_CV20_covar.nlme, main = "PM Sparse (3pt, 20% CV)")
plot(PM_3pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (3pt, 20% CV, weighted)")
Sparse (4pt, 20% CV)
plot(PM_4pt_CV20_covar.nlme, main = "PM Sparse (4pt, 20% CV)")
plot(PM_4pt_CV20_covar.nlme %>% update(weights = varPower()),
main = "PM Sparse (4pt, 20% CV, weighted)")
Adjusting PM Variance Structure
PM_3pt_covar.nlme <- PM_3pt_covar.nlme %>% update(weights = varPower())
PM_4pt_covar.nlme <- PM_4pt_covar.nlme %>% update(weights = varPower())
PM_CV5_covar.nlme <- PM_CV5_covar.nlme %>% update(weights = varPower())
PM_3pt_CV5_covar.nlme <- PM_3pt_CV5_covar.nlme %>% update(weights = varPower())
PM_4pt_CV5_covar.nlme <- PM_4pt_CV5_covar.nlme %>% update(weights = varPower())
PM_CV10_covar.nlme <- PM_CV10_covar.nlme %>% update(weights = varPower())
PM_3pt_CV10_covar.nlme <- PM_3pt_CV10_covar.nlme %>% update(weights = varPower())
PM_4pt_CV10_covar.nlme <- PM_4pt_CV10_covar.nlme %>% update(weights = varPower())
PM_CV20_covar.nlme <- PM_CV20_covar.nlme %>% update(weights = varPower())
PM_3pt_CV20_covar.nlme <- PM_3pt_CV20_covar.nlme %>% update(weights = varPower())
PM_4pt_CV20_covar.nlme <- PM_4pt_CV20_covar.nlme %>% update(weights = varPower())PM Covariate Model Summary Data
# Summary Data Table ------------
PM_covaref_CV.table <- rbind(
fixef(PM_covar.nlme),
fixef(PM_3pt_covar.nlme),
fixef(PM_4pt_covar.nlme),
fixef(PM_CV5_covar.nlme),
fixef(PM_3pt_CV5_covar.nlme),
fixef(PM_4pt_CV5_covar.nlme),
fixef(PM_CV10_covar.nlme),
fixef(PM_3pt_CV10_covar.nlme),
fixef(PM_4pt_CV10_covar.nlme),
fixef(PM_CV20_covar.nlme),
fixef(PM_3pt_CV20_covar.nlme),
fixef(PM_4pt_CV20_covar.nlme)) %>%
as_tibble() %>%
mutate(`Covariate Model` = c("Rich","Sparse (3pt)","Sparse (4pt)",
"Rich (5% CV)","Sparse (5% CV, 3pt)","Sparse (5% CV, 4pt)",
"Rich (10% CV)","Sparse (10% CV, 3pt)","Sparse (10% CV, 4pt)",
"Rich (20% CV)","Sparse (20% CV, 3pt)","Sparse (20% CV, 4pt)")) %>%
pivot_longer(`Vmax.(Intercept)`:`Km.Diplotype4/5`,
names_to = "Parameter",
values_to = "Value") %>%
pivot_wider(names_from = `Covariate Model`, values_from = Value) %>%
separate(Parameter, sep = "\\.", into = c("Parameter","Diplotype"), remove = T) %>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = "Reference (4/41)"),
across(where(is.numeric), round, 2))
# Complete Covariate Effects Table --------------------
PM_covar_Vmax.table <- PM_covaref_CV.table %>%filter(Parameter == "Vmax")
PM_covar_Km.table <- PM_covaref_CV.table %>% filter(Parameter == "Km")
# Vmax Estimates -------------------------------
PM_final_vmax.est <- PM_covar_Vmax.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (4/41)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (4/41)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (4/41)" = "Diplotype4/41")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Km), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Vmax)/Vmax)*100, digits = 2),
Vmax = round(Vmax, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Vmax, Predicted, `Error (%)`) %>%
rename("Predicted (Vmax)" = "Predicted") %>%
arrange(Condition)
# Km Estimates -----------------------------------
PM_final_km.est <- PM_covar_Km.table %>%
pivot_longer(cols = Rich:`Sparse (20% CV, 4pt)`,
names_to = "Condition",
values_to = "Estimate (Eta)") %>%
group_by(Condition) %>%
mutate(Predicted = if_else(Diplotype != "Reference (4/41)",
`Estimate (Eta)` + `Estimate (Eta)`[Diplotype == "Reference (4/41)"],
`Estimate (Eta)`),
Diplotype = recode(Diplotype, "Reference (4/41)" = "Diplotype4/41")) %>%
ungroup() %>%
separate(Diplotype, sep = "Diplotype", remove = T, into = c("Type","Diplotype")) %>%
select(-Type) %>%
left_join(diplotype.data %>% select(-Vmax), by = "Diplotype") %>%
mutate(`Error (%)` = round(abs((Predicted-Km)/Km)*100, digits = 2),
Km = round(Km, digits = 2),
Condition = factor(Condition, levels =
c("Rich",
"Rich (5% CV)",
"Rich (10% CV)",
"Rich (20% CV)",
"Sparse (3pt)",
"Sparse (5% CV, 3pt)",
"Sparse (10% CV, 3pt)",
"Sparse (20% CV, 3pt)",
"Sparse (4pt)",
"Sparse (5% CV, 4pt)",
"Sparse (10% CV, 4pt)",
"Sparse (20% CV, 4pt)"))) %>%
select(Parameter, Diplotype, Condition, Km, Predicted, `Error (%)`) %>%
rename("Predicted (Km)" = "Predicted") %>%
arrange(Condition)Model Accuracy
Combining Data Across All Study Design
all_km.est <- rbind(final_km.est, PM_final_km.est)
all_vmax.est <- rbind(final_vmax.est, PM_final_vmax.est)Predicted vs Actual Figure (Vmax)
Vmax.plot <- ggplot(data = all_vmax.est, aes(x = Vmax, y = `Predicted (Vmax)`))+
geom_point(size = 3, alpha = 0.85, aes(fill = Diplotype), shape = 21)+
geom_abline(slope = 1, intercept = 0, color = "red")+
geom_abline(slope = 1, intercept = 0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.2, color = "red", linetype = "dashed")+
scale_y_log10()+
scale_x_log10()+
facet_wrap(.~Condition, ncol = 4, scales = "free")+
xlab("Actual (Vmax)")+
theme_bw(base_size = 12)+
ggeasy::easy_move_legend("top")+
scale_fill_viridis_d()
Vmax.plot
Predicted vs Actual Figure (Km)
Km.plot <- ggplot(data = all_km.est, aes(x = Km, y = `Predicted (Km)`))+
geom_point(size = 3, alpha = 0.85, aes(fill = Diplotype), shape = 21)+
geom_abline(slope = 1, intercept = 0, color = "red")+
geom_abline(slope = 1, intercept = 0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.2, color = "red", linetype = "dashed")+
scale_y_log10()+
scale_x_log10()+
facet_wrap(.~Condition, ncol = 4, scales = "free")+
xlab("Actual (Km)")+
theme_bw(base_size = 12)+
ggeasy::easy_move_legend("top")+
scale_fill_viridis_d()
Km.plot
Summary Tables
t-Table Extractor Function
# Function to grab t-Table output from NLME model objects -----------
get_tTable <- function(data.nlme, condition, PM = FALSE){
label <- if_else(PM == TRUE, "4/41", "1/1")
summary(data.nlme)$tTable %>%
as_tibble() %>%
mutate(Parameter = row.names(summary(data.nlme)$tTable)) %>%
relocate(Parameter) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = paste0("Diplotype",label))) %>%
group_by(Parameter) %>%
mutate(Value = if_else(Diplotype == paste0("Diplotype",label),
Value,
Value+Value[Diplotype == paste0("Diplotype",label)])) %>%
separate(Diplotype, remove = T, sep = "Diplotype", into = c(".", "Diplotype")) %>%
select(Parameter, Diplotype, Value, Std.Error) %>%
mutate(`RSE(%)` = round((Std.Error/Value)*100, digits = 2),
Condition = condition)
}Final Model Parameter Estimates and Standard Error Datasets
# Standard Error Tables
## EM and UM Population -------
rich_CV0_covar.tTable <- get_tTable(rich_covar.nlme, "Rich (0% CV)")
sparse3pt_CV0_covar.tTable <- get_tTable(sparse3pt_covar.nlme, "Sparse (0% CV, 3pt)")
sparse4pt_CV0_covar.tTable <- get_tTable(sparse4pt_covar.nlme, "Sparse (0% CV, 4pt)")
rich_CV5_covar.tTable <- get_tTable(rich_CV5_covar.nlme, "Rich (5% CV)")
sparse3pt_CV5_covar.tTable <- get_tTable(sparse3pt_CV5_covar.nlme, "Sparse (5% CV, 3pt)")
sparse4pt_CV5_covar.tTable <- get_tTable(sparse4pt_CV5_covar.nlme, "Sparse (5% CV, 4pt)")
rich_CV10_covar.tTable <- get_tTable(rich_CV10_covar.nlme, "Rich (10% CV)")
sparse3pt_CV10_covar.tTable <- get_tTable(sparse3pt_CV10_covar.nlme, "Sparse (10% CV, 3pt)")
sparse4pt_CV10_covar.tTable <- get_tTable(sparse4pt_CV10_covar.nlme, "Sparse (10% CV, 4pt)")
rich_CV20_covar.tTable <- get_tTable(rich_CV20_covar.nlme, "Rich (20% CV)")
sparse3pt_CV20_covar.tTable <- get_tTable(sparse3pt_CV20_covar.nlme, "Sparse (20% CV, 3pt)")
sparse4pt_CV20_covar.tTable <- get_tTable(sparse4pt_CV20_covar.nlme, "Sparse (20% CV, 4pt)")
## PM Population -------
PM_CV0_covar.tTable <- get_tTable(PM_covar.nlme, "Rich (0% CV)", T)
PM_3pt_CV0_covar.tTable <- get_tTable(PM_3pt_covar.nlme, "Sparse (0% CV, 3pt)", T)
PM_4pt_CV0_covar.tTable <- get_tTable(PM_4pt_covar.nlme, "Sparse (0% CV, 4pt)", T)
PM_CV5_covar.tTable <- get_tTable(PM_CV5_covar.nlme, "Rich (5% CV)", T)
PM_3pt_CV5_covar.tTable <- get_tTable(PM_3pt_CV5_covar.nlme, "Sparse (5% CV, 3pt)", T)
PM_4pt_CV5_covar.tTable <- get_tTable(PM_4pt_CV5_covar.nlme, "Sparse (5% CV, 4pt)", T)
PM_CV10_covar.tTable <- get_tTable(PM_CV10_covar.nlme, "Rich (10% CV)", T)
PM_3pt_CV10_covar.tTable <- get_tTable(PM_3pt_CV10_covar.nlme, "Sparse (10% CV, 3pt)", T)
PM_4pt_CV10_covar.tTable <- get_tTable(PM_4pt_CV10_covar.nlme, "Sparse (10% CV, 4pt)", T)
PM_CV20_covar.tTable <- get_tTable(PM_CV20_covar.nlme, "Rich (20% CV)", T)
PM_3pt_CV20_covar.tTable <- get_tTable(PM_3pt_CV20_covar.nlme, "Sparse (20% CV, 3pt)", T)
PM_4pt_CV20_covar.tTable <- get_tTable(PM_4pt_CV20_covar.nlme, "Sparse (20% CV, 4pt)", T)Combining Datasets across Experimental Conditions
## Combined tTable Outpout -----------------
complete_CV_covar.tTable <- rbind(
# EM and UM Population Estimates -------
rich_CV5_covar.tTable,
sparse3pt_CV5_covar.tTable,
sparse4pt_CV5_covar.tTable,
rich_CV10_covar.tTable,
sparse3pt_CV10_covar.tTable,
sparse4pt_CV10_covar.tTable,
rich_CV20_covar.tTable,
sparse3pt_CV20_covar.tTable,
sparse4pt_CV20_covar.tTable,
#PM and IM Population Estimates -----
PM_CV5_covar.tTable,
PM_3pt_CV5_covar.tTable,
PM_4pt_CV5_covar.tTable,
PM_CV10_covar.tTable,
PM_3pt_CV10_covar.tTable,
PM_4pt_CV10_covar.tTable,
PM_CV20_covar.tTable,
PM_3pt_CV20_covar.tTable,
PM_4pt_CV20_covar.tTable)Confidence Interval Data (All Conditions)
# Confidence Interval Table ---------
## EM and UM Population Intervals
rich_CV5_covar.intervals <- intervals(rich_CV5_covar.nlme, which = "fixed")
sparse3pt_CV5_covar.intervals <- intervals(sparse3pt_CV5_covar.nlme, which = "fixed")
sparse4pt_CV5_covar.intervals <- intervals(sparse4pt_CV5_covar.nlme, which = "fixed")
rich_CV10_covar.intervals <- intervals(rich_CV10_covar.nlme, which = "fixed")
sparse3pt_CV10_covar.intervals <- intervals(sparse3pt_CV10_covar.nlme, which = "fixed")
sparse4pt_CV10_covar.intervals <- intervals(sparse4pt_CV10_covar.nlme, which = "fixed")
rich_CV20_covar.intervals <- intervals(rich_CV20_covar.nlme, which = "fixed")
sparse3pt_CV20_covar.intervals <- intervals(sparse3pt_CV20_covar.nlme, which = "fixed")
sparse4pt_CV20_covar.intervals <- intervals(sparse4pt_CV20_covar.nlme, which = "fixed")
## PM and IM Population Intervals
PM_covar.intervals <- intervals(PM_covar.nlme, which = "fixed")
PM_3pt_covar.intervals <- intervals(PM_3pt_covar.nlme, which = "fixed")
PM_4pt_covar.intervals <- intervals(PM_4pt_covar.nlme, which = "fixed")
PM_CV5_covar.intervals <- intervals(PM_CV5_covar.nlme, which = "fixed")
PM_3pt_CV5_covar.intervals <- intervals(PM_3pt_CV5_covar.nlme, which = "fixed")
PM_4pt_CV5_covar.intervals <- intervals(PM_4pt_CV5_covar.nlme, which = "fixed")
PM_CV10_covar.intervals <- intervals(PM_CV10_covar.nlme, which = "fixed")
PM_3pt_CV10_covar.intervals <- intervals(PM_3pt_CV10_covar.nlme, which = "fixed")
PM_4pt_CV10_covar.intervals <- intervals(PM_4pt_CV10_covar.nlme, which = "fixed")
PM_CV20_covar.intervals <- intervals(PM_CV20_covar.nlme, which = "fixed")
PM_3pt_CV20_covar.intervals <- intervals(PM_3pt_CV20_covar.nlme, which = "fixed")
PM_4pt_CV20_covar.intervals <- intervals(PM_4pt_CV20_covar.nlme, which = "fixed")Confidence Interval Extractor Function
# Function to extract and tidy confidence interval estimates ---------
extract_CI3 <- function(int.data, name, PM = FALSE){
label <- if_else(PM == TRUE, "4/41", "1/1")
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(Parameter = rownames(int.data$fixed)) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = paste0("Diplotype",label))) %>%
group_by(Parameter) %>%
mutate(lower = if_else(Diplotype != paste0("Diplotype",label),
lower + est.[Diplotype == paste0("Diplotype",label)],
lower),
upper = if_else(Diplotype != paste0("Diplotype",label),
upper + est.[Diplotype == paste0("Diplotype",label)],
upper),
est. = if_else(Diplotype != paste0("Diplotype",label),
est. + est.[Diplotype == paste0("Diplotype",label)],
est.),
across(where(is.numeric), round, 2),
`CI (95%)` = paste0("(",lower,", ",upper,")")) %>%
ungroup() %>%
separate(Diplotype, remove = T, sep = "Diplotype", into = c(".", "Diplotype"))%>%
select(Parameter, Diplotype, est., `CI (95%)`)
colnames(temp) <- c("Parameter", "Diplotype",
paste(name, " Estimate"),
paste(name, " CI (95%)"))
temp
}Generating Confidence Interval Datatables
# Combined Tables
CV0_covar_ci.table <- cbind(
extract_CI3(rich_covar.intervals, name = "Rich (0% CV)"),
rich_CV0_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_covar.intervals, name = "Sparse 3pt (0% CV)")[,3:4],
sparse3pt_CV0_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_covar.intervals, name = "Sparse 4pt (0% CV)")[,3:4],
sparse4pt_CV0_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
CV5_covar_ci.table <- cbind(
extract_CI3(rich_CV5_covar.intervals, name = "Rich (5% CV)"),
rich_CV5_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_CV5_covar.intervals, name = "Sparse 3pt (5% CV)")[,3:4],
sparse3pt_CV5_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_CV5_covar.intervals, name = "Sparse 4pt (5% CV)")[,3:4],
sparse4pt_CV5_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
CV10_covar_ci.table <- cbind(
extract_CI3(rich_CV10_covar.intervals, name = "Rich (10% CV)"),
rich_CV10_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_CV10_covar.intervals, name = "Sparse 3pt (10% CV)")[,3:4],
sparse3pt_CV10_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_CV10_covar.intervals, name = "Sparse 4pt (10% CV)")[,3:4],
sparse4pt_CV10_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
CV20_covar_ci.table <- cbind(
extract_CI3(rich_CV20_covar.intervals, name = "Rich (20% CV)"),
rich_CV20_covar.tTable$`RSE(%)`,
extract_CI3(sparse3pt_CV20_covar.intervals, name = "Sparse 3pt (20% CV)")[,3:4],
sparse3pt_CV20_covar.tTable$`RSE(%)`,
extract_CI3(sparse4pt_CV20_covar.intervals, name = "Sparse 4pt (20% CV)")[,3:4],
sparse4pt_CV20_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV0_covar_ci.table <- cbind(
extract_CI3(PM_covar.intervals, name = "Rich (0% CV)", T),
PM_CV0_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_covar.intervals, name = "Sparse 3pt (0% CV)", T)[,3:4],
PM_3pt_CV0_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_covar.intervals, name = "Sparse 4pt (0% CV)", T)[,3:4],
PM_4pt_CV0_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV5_covar_ci.table <- cbind(
extract_CI3(PM_CV5_covar.intervals, name = "Rich (5% CV)", T),
PM_CV5_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_CV5_covar.intervals, name = "Sparse 3pt (5% CV)", T)[,3:4],
PM_3pt_CV5_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_CV5_covar.intervals, name = "Sparse 4pt (5% CV)", T)[,3:4],
PM_4pt_CV5_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV10_covar_ci.table <- cbind(
extract_CI3(PM_CV10_covar.intervals, name = "Rich (10% CV)", T),
PM_CV10_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_CV10_covar.intervals, name = "Sparse 3pt (10% CV)", T)[,3:4],
PM_3pt_CV10_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_CV10_covar.intervals, name = "Sparse 4pt (10% CV)", T)[,3:4],
PM_4pt_CV10_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)
PM_CV20_covar_ci.table <- cbind(
extract_CI3(PM_CV20_covar.intervals, name = "Rich (20% CV)", T),
PM_CV20_covar.tTable$`RSE(%)`,
extract_CI3(PM_3pt_CV20_covar.intervals, name = "Sparse 3pt (20% CV)", T)[,3:4],
PM_3pt_CV20_covar.tTable$`RSE(%)`,
extract_CI3(PM_4pt_CV20_covar.intervals, name = "Sparse 4pt (20% CV)", T)[,3:4],
PM_4pt_CV20_covar.tTable$`RSE(%)`) %>%
left_join(actual.est, by = c("Diplotype","Parameter")) %>%
relocate(Actual, .after = Diplotype)Renaming Datatable Columns
colnames(CV0_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(CV5_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(CV10_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(CV20_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV0_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV5_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV10_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")
colnames(PM_CV20_covar_ci.table) <- c("Parameter","Diplotype","Actual",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)",
"Estimate", "CI (95%)", "RSE(%)")Covariate Model CI Tables (UM/EM)
Covariate Model (0% CV) CI Table
CV0_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 0% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich" = 3,
"Sparse (3pt)" = 3,
"Sparse (4pt)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich
|
Sparse (3pt)
|
Sparse (4pt)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 367.01 | (319.69, 414.33) | 6.64 | 367.01 | (319.34, 414.67) | 6.79 | 367.01 | (319.51, 414.51) | 6.73 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 486.93 | (420.01, 553.85) | 7.08 | 486.93 | (419.52, 554.34) | 7.24 | 486.93 | (419.75, 554.11) | 7.18 |
| Vmax | 1/2x2 | 1410.00 | 1459.75 | (1392.84, 1526.67) | 2.36 | 1459.75 | (1392.34, 1527.16) | 2.42 | 1459.75 | (1392.57, 1526.93) | 2.39 |
| Vmax | 1/4 | 274.00 | 283.67 | (216.75, 350.58) | 12.14 | 283.67 | (216.26, 351.08) | 12.43 | 283.67 | (216.49, 350.85) | 12.32 |
| Vmax | 2/2 | 252.00 | 260.89 | (193.98, 327.81) | 13.21 | 260.89 | (193.48, 328.3) | 13.52 | 260.89 | (193.71, 328.07) | 13.40 |
| Vmax | 2/3 | 251.00 | 259.86 | (192.94, 326.77) | 13.26 | 259.86 | (192.45, 327.27) | 13.57 | 259.86 | (192.68, 327.04) | 13.45 |
| Vmax | 2/4 | 95.90 | 99.28 | (32.37, 166.2) | 34.70 | 99.28 | (31.87, 166.69) | 35.52 | 99.28 | (32.1, 166.46) | 35.20 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.50 | (1.33, 1.67) | 5.92 | 1.50 | (1.33, 1.67) | 6.06 | 1.50 | (1.33, 1.67) | 6.00 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.44 | (1.2, 1.68) | 8.73 | 1.44 | (1.19, 1.68) | 8.93 | 1.44 | (1.19, 1.68) | 8.85 |
| Km | 1/2x2 | 1.82 | 1.89 | (1.65, 2.13) | 6.65 | 1.89 | (1.64, 2.14) | 6.80 | 1.89 | (1.64, 2.13) | 6.74 |
| Km | 1/4 | 0.92 | 0.95 | (0.71, 1.2) | 13.18 | 0.95 | (0.71, 1.2) | 13.49 | 0.95 | (0.71, 1.2) | 13.37 |
| Km | 2/2 | 4.59 | 4.75 | (4.51, 5) | 2.64 | 4.75 | (4.51, 5) | 2.71 | 4.75 | (4.51, 5) | 2.68 |
| Km | 2/3 | 1.32 | 1.37 | (1.12, 1.61) | 9.19 | 1.37 | (1.12, 1.61) | 9.41 | 1.37 | (1.12, 1.61) | 9.32 |
| Km | 2/4 | 1.68 | 1.74 | (1.49, 1.98) | 7.24 | 1.74 | (1.49, 1.98) | 7.41 | 1.74 | (1.49, 1.98) | 7.34 |
Covariate Model (5% CV) CI Table
CV5_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 5% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (5% CV)" = 3,
"Sparse (3pt, 5% CV)" = 3,
"Sparse (4pt, 5% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich (5% CV)
|
Sparse (3pt, 5% CV)
|
Sparse (4pt, 5% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 362.01 | (326.7, 397.31) | 5.02 | 364.79 | (346.49, 383.1) | 2.63 | 363.04 | (343.06, 383.02) | 2.86 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 487.84 | (437.63, 538.06) | 5.30 | 482.00 | (452.75, 511.25) | 3.17 | 481.46 | (451.75, 511.16) | 3.21 |
| Vmax | 1/2x2 | 1410.00 | 1461.41 | (1405.98, 1516.83) | 1.95 | 1463.40 | (1387.76, 1539.04) | 2.70 | 1460.85 | (1405.71, 1515.98) | 1.96 |
| Vmax | 1/4 | 274.00 | 283.50 | (233.72, 333.28) | 9.04 | 278.28 | (254.56, 302.01) | 4.46 | 278.34 | (250.88, 305.8) | 5.13 |
| Vmax | 2/2 | 252.00 | 261.57 | (211.72, 311.41) | 9.81 | 264.06 | (239.64, 288.48) | 4.84 | 264.62 | (236.97, 292.26) | 5.44 |
| Vmax | 2/3 | 251.00 | 261.41 | (211.64, 311.17) | 9.80 | 258.97 | (235.43, 282.51) | 4.75 | 259.09 | (231.73, 286.45) | 5.49 |
| Vmax | 2/4 | 95.90 | 97.83 | (48.22, 147.45) | 26.11 | 98.79 | (77.14, 120.44) | 11.46 | 98.20 | (71.53, 124.86) | 14.13 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.47 | (1.34, 1.59) | 4.26 | 1.49 | (1.38, 1.59) | 3.74 | 1.47 | (1.36, 1.58) | 3.87 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.46 | (1.29, 1.63) | 6.03 | 1.46 | (1.31, 1.62) | 5.41 | 1.47 | (1.31, 1.62) | 5.54 |
| Km | 1/2x2 | 1.82 | 1.88 | (1.71, 2.06) | 4.81 | 1.89 | (1.71, 2.08) | 5.18 | 1.90 | (1.72, 2.07) | 4.87 |
| Km | 1/4 | 0.92 | 0.97 | (0.8, 1.13) | 8.93 | 0.96 | (0.83, 1.09) | 6.97 | 0.94 | (0.79, 1.08) | 7.96 |
| Km | 2/2 | 4.59 | 4.83 | (4.61, 5.05) | 2.37 | 4.99 | (4.65, 5.32) | 3.51 | 4.99 | (4.73, 5.25) | 2.71 |
| Km | 2/3 | 1.32 | 1.39 | (1.22, 1.56) | 6.33 | 1.37 | (1.22, 1.51) | 5.48 | 1.37 | (1.22, 1.52) | 5.72 |
| Km | 2/4 | 1.68 | 1.68 | (1.5, 1.85) | 5.31 | 1.67 | (1.52, 1.81) | 4.61 | 1.65 | (1.5, 1.8) | 4.73 |
Covariate Model (10% CV) CI Table
CV10_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 10% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (10% CV)" = 3,
"Sparse (3pt, 10% CV)" = 3,
"Sparse (4pt, 10% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich (10% CV)
|
Sparse (3pt, 10% CV)
|
Sparse (4pt, 10% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 355.96 | (336.28, 375.64) | 2.85 | 365.54 | (341.6, 389.48) | 3.43 | 361.21 | (339.36, 383.05) | 3.15 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 487.46 | (457.38, 517.54) | 3.18 | 481.34 | (442.9, 519.79) | 4.18 | 479.61 | (445.64, 513.58) | 3.68 |
| Vmax | 1/2x2 | 1410.00 | 1479.15 | (1420.74, 1537.56) | 2.03 | 1476.20 | (1380.61, 1571.8) | 3.39 | 1474.03 | (1400.08, 1547.98) | 2.61 |
| Vmax | 1/4 | 274.00 | 282.04 | (255.33, 308.75) | 4.88 | 273.99 | (243.4, 304.57) | 5.84 | 275.03 | (246.08, 303.98) | 5.48 |
| Vmax | 2/2 | 252.00 | 263.89 | (236.56, 291.22) | 5.33 | 267.54 | (235.45, 299.63) | 6.27 | 268.41 | (238.59, 298.23) | 5.78 |
| Vmax | 2/3 | 251.00 | 262.09 | (235.47, 288.71) | 5.23 | 259.86 | (229.39, 290.33) | 6.13 | 259.08 | (230.28, 287.89) | 5.78 |
| Vmax | 2/4 | 95.90 | 96.22 | (70.87, 121.57) | 13.57 | 98.98 | (71.55, 126.42) | 14.50 | 97.58 | (70.7, 124.47) | 14.33 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.43 | (1.33, 1.53) | 3.63 | 1.48 | (1.32, 1.65) | 5.85 | 1.45 | (1.3, 1.61) | 5.43 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.49 | (1.34, 1.63) | 5.03 | 1.51 | (1.28, 1.75) | 8.22 | 1.51 | (1.29, 1.72) | 7.54 |
| Km | 1/2x2 | 1.82 | 1.91 | (1.74, 2.07) | 4.51 | 1.92 | (1.64, 2.2) | 7.51 | 1.92 | (1.67, 2.17) | 6.74 |
| Km | 1/4 | 0.92 | 0.97 | (0.85, 1.1) | 6.55 | 0.96 | (0.76, 1.17) | 10.98 | 0.96 | (0.78, 1.15) | 10.12 |
| Km | 2/2 | 4.59 | 4.97 | (4.63, 5.31) | 3.54 | 5.19 | (4.64, 5.74) | 5.54 | 5.19 | (4.7, 5.69) | 4.98 |
| Km | 2/3 | 1.32 | 1.41 | (1.26, 1.55) | 5.17 | 1.38 | (1.15, 1.6) | 8.61 | 1.37 | (1.16, 1.57) | 7.95 |
| Km | 2/4 | 1.68 | 1.61 | (1.46, 1.76) | 4.76 | 1.62 | (1.38, 1.86) | 7.70 | 1.59 | (1.37, 1.8) | 7.19 |
Covariate Model (20% CV) CI Table
CV20_covar_ci.table %>%
kbl(caption = "EM and UM Covariate Model Estimates (w/ 20% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (20% CV)" = 3,
"Sparse (3pt, 20% CV)" = 3,
"Sparse (4pt, 20% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Wild-Type)", 1,1) %>%
pack_rows("Variant Estimates", 2, 7)%>%
pack_rows("Km Estimate (Wild-Type)", 8,8) %>%
pack_rows("Variant Estimates", 9, 14)|
Rich (20% CV)
|
Sparse (3pt, 20% CV)
|
Sparse (4pt, 20% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Wild-Type) | |||||||||||
| Vmax | 1/1 | 354.50 | 345.88 | (319.99, 371.77) | 3.85 | 366.90 | (329.06, 404.73) | 5.39 | 358.47 | (326.75, 390.19) | 4.60 |
| Variant Estimates | |||||||||||
| Vmax | 1/2 | 470.33 | 491.31 | (448.44, 534.19) | 4.49 | 484.18 | (421.45, 546.91) | 6.78 | 477.51 | (426, 529.02) | 5.61 |
| Vmax | 1/2x2 | 1410.00 | 1510.05 | (1408.03, 1612.07) | 3.48 | 1501.07 | (1339.29, 1662.86) | 5.64 | 1497.51 | (1373.43, 1621.59) | 4.31 |
| Vmax | 1/4 | 274.00 | 280.85 | (246.74, 314.97) | 6.25 | 265.12 | (218.62, 311.62) | 9.18 | 267.56 | (227.18, 307.95) | 7.85 |
| Vmax | 2/2 | 252.00 | 266.92 | (230.87, 302.98) | 6.95 | 275.51 | (224.49, 326.53) | 9.69 | 275.49 | (232.06, 318.91) | 8.20 |
| Vmax | 2/3 | 251.00 | 264.38 | (230.39, 298.38) | 6.62 | 262.08 | (215.29, 308.86) | 9.34 | 259.33 | (218.97, 299.7) | 8.10 |
| Vmax | 2/4 | 95.90 | 93.65 | (63.33, 123.97) | 16.67 | 98.20 | (58.18, 138.22) | 21.32 | 95.77 | (60.19, 131.34) | 19.32 |
| Km Estimate (Wild-Type) | |||||||||||
| Km | 1/1 | 1.45 | 1.37 | (1.18, 1.55) | 6.92 | 1.48 | (1.15, 1.8) | 11.68 | 1.43 | (1.14, 1.72) | 10.53 |
| Variant Estimates | |||||||||||
| Km | 1/2 | 1.39 | 1.56 | (1.28, 1.83) | 9.13 | 1.64 | (1.16, 2.13) | 15.47 | 1.61 | (1.18, 2.04) | 13.84 |
| Km | 1/2x2 | 1.82 | 1.94 | (1.63, 2.25) | 8.25 | 1.98 | (1.45, 2.51) | 14.06 | 1.98 | (1.51, 2.45) | 12.32 |
| Km | 1/4 | 0.92 | 1.00 | (0.77, 1.23) | 11.91 | 0.96 | (0.54, 1.37) | 22.61 | 0.98 | (0.62, 1.34) | 19.16 |
| Km | 2/2 | 4.59 | 5.20 | (4.52, 5.88) | 6.71 | 5.62 | (4.48, 6.76) | 10.62 | 5.58 | (4.53, 6.64) | 9.82 |
| Km | 2/3 | 1.32 | 1.45 | (1.18, 1.72) | 9.47 | 1.41 | (0.95, 1.86) | 17.00 | 1.38 | (0.98, 1.78) | 15.24 |
| Km | 2/4 | 1.68 | 1.51 | (1.24, 1.79) | 9.29 | 1.53 | (1.06, 2) | 16.19 | 1.46 | (1.04, 1.88) | 14.96 |
Covariate Model CI Tables (IM/PM)
IM/PM Covariate Model (0% CV) CI Table
PM_CV0_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 0% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich" = 3,
"Sparse (3pt)" = 3,
"Sparse (4pt)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich
|
Sparse (3pt)
|
Sparse (4pt)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 31.99 | (30.14, 33.84) | 2.96 | 31.99 | (30.28, 33.69) | 2.72 | 31.94 | (30.2, 33.67) | 2.79 |
| Vmax | 4/5 | 12.32 | 12.75 | (10.14, 15.37) | 10.48 | 12.74 | (10.34, 15.15) | 9.65 | 12.43 | (10.64, 14.23) | 7.40 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.54 | (1.69, 9.39) | 35.57 | 5.61 | (5.3, 5.92) | 2.83 | 5.54 | (5.07, 6.02) | 4.38 |
| Km | 4/5 | 69.15 | 71.60 | (66.15, 77.05) | 3.89 | 71.96 | (67.76, 76.15) | 2.98 | 68.35 | (65.66, 71.04) | 2.02 |
IM/PM Covariate Model (5% CV) CI Table
PM_CV5_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 5% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (5% CV)" = 3,
"Sparse (3pt, 5% CV)" = 3,
"Sparse (4pt, 5% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich (5% CV)
|
Sparse (3pt, 5% CV)
|
Sparse (4pt, 5% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 31.70 | (30.04, 33.35) | 2.67 | 31.88 | (30.1, 33.65) | 2.85 | 31.62 | (29.89, 33.36) | 2.81 |
| Vmax | 4/5 | 12.32 | 12.60 | (10.27, 14.92) | 9.44 | 12.47 | (10.11, 14.84) | 9.70 | 12.56 | (10.72, 14.4) | 7.51 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.56 | (5.17, 5.96) | 3.65 | 6.19 | (5.07, 7.31) | 9.23 | 5.57 | (5.03, 6.11) | 4.97 |
| Km | 4/5 | 69.15 | 69.17 | (67.09, 71.25) | 1.54 | 67.23 | (61.18, 73.28) | 4.60 | 69.26 | (65.67, 72.85) | 2.66 |
PM Covariate Model (10% CV) CI Table
PM_CV10_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 10% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (10% CV)" = 3,
"Sparse (3pt, 10% CV)" = 3,
"Sparse (4pt, 10% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich (10% CV)
|
Sparse (3pt, 10% CV)
|
Sparse (4pt, 10% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 31.37 | (29.86, 32.87) | 2.45 | 32.02 | (29.53, 34.52) | 3.98 | 31.40 | (29.23, 33.57) | 3.55 |
| Vmax | 4/5 | 12.32 | 12.50 | (10.44, 14.56) | 8.44 | 12.25 | (9.58, 14.92) | 11.12 | 12.46 | (10.1, 14.81) | 9.68 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.60 | (5.2, 5.99) | 3.60 | 6.82 | (4.18, 9.46) | 19.74 | 5.66 | (4.87, 6.45) | 7.16 |
| Km | 4/5 | 69.15 | 67.53 | (64.29, 70.77) | 2.46 | 63.57 | (54.88, 72.26) | 6.99 | 68.20 | (62.5, 73.89) | 4.28 |
PM Covariate Model (20% CV) CI Table
PM_CV20_covar_ci.table %>%
kbl(caption = "IM and PM Covariate Model Estimates (w/ 20% Residual Error & 95% CI)") %>%
kable_classic(full_width = T, "striped")%>%
add_header_above(c(" " = 1," " = 1, " " = 1,
"Rich (20% CV)" = 3,
"Sparse (3pt, 20% CV)" = 3,
"Sparse (4pt, 20% CV)" = 3)) %>%
pack_rows("Vmax Estimate (Variant)", 1,2) %>%
pack_rows("Km Estimate (Variant)", 3,4)|
Rich (20% CV)
|
Sparse (3pt, 20% CV)
|
Sparse (4pt, 20% CV)
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Parameter | Diplotype | Actual | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) | Estimate | CI (95%) | RSE(%) |
| Vmax Estimate (Variant) | |||||||||||
| Vmax | 4/41 | 30.90 | 30.77 | (29.18, 32.36) | 2.64 | 32.17 | (28.23, 36.11) | 6.26 | 31.02 | (27.58, 34.46) | 5.68 |
| Vmax | 4/5 | 12.32 | 12.33 | (10.37, 14.3) | 8.16 | 11.89 | (7.68, 16.1) | 18.10 | 12.28 | (8.52, 16.03) | 15.67 |
| Km Estimate (Variant) | |||||||||||
| Km | 4/41 | 5.35 | 5.70 | (4.93, 6.48) | 6.96 | 8.33 | (3.77, 12.89) | 27.95 | 5.85 | (4.43, 7.28) | 12.49 |
| Km | 4/5 | 69.15 | 64.82 | (59.1, 70.55) | 4.51 | 57.03 | (42.81, 71.26) | 12.74 | 65.65 | (55.21, 76.09) | 8.15 |
Final Model Predictions
NLME Predictor Function
The predictor function sim_nlme() was scripted to
generate tidy data frames of NLME predictions given the specified
substrate concentration range. The inputs for this function are as
followed:
tTable - extracted NLME model tTable estimates (see Final Model Parameter Estimates and Standard Error Datasets section).
label - Specifies the condition being simulated (ex. “Rich 0% CV”).
Group - Specifies the CYP2D6 metabolizer category (ex. “IM & PM”).
Start - Substrate concentration to start simulation (default = 0 uM).
End - Last substrate concentration (default = 100 uM).
by - Simulated concentration increment (default = 0.1)
sim_nlme <- function(tTable, label, Group, Start = 0, End = 100, by = 0.1){
diplotypes_contained <- unique(tTable$Diplotype)
# Simulate Diplotype Rates -----
tTable %>%
select(Parameter, Diplotype, Value) %>%
pivot_wider(names_from = Parameter,
values_from = Value) %>%
expand_grid(S = seq(Start,End,by)) %>%
mutate(V = (Vmax*S)/(Km+S),
Condition = label,
Group = Group)
}Update Data Function (Extended)
update_data_full() is an extended version of the
previous update_data() function; it differs by retaining Km
and Vmax information in the output. It was created as an intermediate
function for the revised update_data2() function (see
below).
update_data_full <-function(`CV%`, type, n_indiv = 10, `CV%.D` = 25, by = 0.1,
start = 0, end = 100, Seed = 23457, Pop.freq = FALSE){
# Sampling Information --------------
# Full Range Set
full_set <- c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)
PM_range <- c(1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, 2000)
# Strategic Sampling Sets ----
sparse3pt_set1 <- c(0.1, 2.5, 50)
sparse3pt_set2 <- c(1, 10, 100)
sparse4pt_set1 <- c(0.1, 2.5, 25, 50)
sparse4pt_set2 <- c(1, 10, 50, 100)
PM_3pt_set1 <- c(10,100,1000)
PM_3pt_set2 <- c(25,250,2000)
PM_4pt_set1 <- c(1,10,100,1000)
PM_4pt_set2 <- c(5,25,250,2000)
PM_range_set1 <- c(1,10,25,50,100,400,1000,2000)
PM_range_set2 <- c(5,15,30,60,90,200,800,1600)
# Rich Sampling Simulation ----------------
Rich <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
# Sparse Sampling Simulation (3 point) ------------
Sparse3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse3pt_set1, S %in% sparse3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
# Sparse Sampling Simulation (4 point) --------------
Sparse4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = end, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% full_set) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% sparse4pt_set1, S %in% sparse4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
PM <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
PM_3pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_3pt_set1, S %in% PM_3pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
PM_4pt <- population.sim(n = n_indiv, `CV%_D` = `CV%.D`, `CV%_V` = `CV%`, By = by,
Start = start, End = 2000, seed = Seed, pop_freq = Pop.freq) %>%
filter(S %in% c(PM_range)) %>%
mutate(Set = if_else(as.numeric(ID) %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1", S %in% PM_4pt_set1, S %in% PM_4pt_set2)) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, Vmax, Km, V, S)
# Output Selection ----------
switch(type,
"rich" = Rich,
"sparse3pt" = Sparse3pt,
"sparse4pt" = Sparse4pt,
"PM" = PM,
"PM_3pt" = PM_3pt,
"PM_4pt" = PM_4pt)
}Update Data Function (Version 2)
Revised vs of update_data() designed to allow user to
include more information (i.e., Condition, and Group) to the simulated
data output.
update_data2 <- function(`CV%`, type, Condition, Group,
n_indiv = 10,
`CV%.D` = 25,
by = 0.1,
start = 0,
end = 100,
Seed = 23457,
Pop.freq = FALSE){
update_data_full(`CV%`, type, n_indiv = n_indiv,
`CV%.D` = `CV%.D`, by = by,
start = start, end = end,
Seed = Seed, Pop.freq = Pop.freq) %>%
mutate(Condition = Condition,
Group = Group) %>%
filter(if_else(Group == "EM & UM",
!Diplotype %in% c("4/41", "4/5"),
Diplotype %in% c("4/41", "4/5")))
}Sensitivity Analysis Data (Population Level)
# Population Level Simulated Data Based on Experimental Conditions ---
sensitivity.data <- rbind(sim_nlme(rich_CV0_covar.tTable,"Rich (0% CV)", "EM & UM"),
sim_nlme(rich_CV5_covar.tTable,"Rich (5% CV)", "EM & UM"),
sim_nlme(rich_CV10_covar.tTable,"Rich (10% CV)", "EM & UM"),
sim_nlme(rich_CV20_covar.tTable,"Rich (20% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV0_covar.tTable,"Sparse (3pt, 0% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV5_covar.tTable,"Sparse (3pt, 5% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV10_covar.tTable,"Sparse (3pt, 10% CV)", "EM & UM"),
sim_nlme(sparse3pt_CV20_covar.tTable,"Sparse (3pt, 20% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV0_covar.tTable,"Sparse (4pt, 0% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV5_covar.tTable,"Sparse (4pt, 5% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV10_covar.tTable,"Sparse (4pt, 10% CV)", "EM & UM"),
sim_nlme(sparse4pt_CV20_covar.tTable,"Sparse (4pt, 20% CV)", "EM & UM"),
sim_nlme(PM_CV0_covar.tTable, "Rich (0% CV)", "IM & PM", End = 2000),
sim_nlme(PM_CV5_covar.tTable, "Rich (5% CV)", "IM & PM", End = 2000),
sim_nlme(PM_CV10_covar.tTable, "Rich (10% CV)", "IM & PM", End = 2000),
sim_nlme(PM_CV20_covar.tTable, "Rich (20% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV0_covar.tTable, "Sparse (3pt, 0% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV5_covar.tTable, "Sparse (3pt, 5% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV10_covar.tTable, "Sparse (3pt, 10% CV)", "IM & PM", End = 2000),
sim_nlme(PM_3pt_CV20_covar.tTable, "Sparse (3pt, 20% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV0_covar.tTable, "Sparse (4pt, 0% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV5_covar.tTable, "Sparse (4pt, 5% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV10_covar.tTable, "Sparse (4pt, 10% CV)", "IM & PM", End = 2000),
sim_nlme(PM_4pt_CV20_covar.tTable, "Sparse (4pt, 20% CV)", "IM & PM", End = 2000))
sensitivity.data$Condition <- factor(sensitivity.data$Condition,
levels = unique(sensitivity.data$Condition))Sensitivity Analysis Data (Individual Level)
# Individual Level Simulated Data Based on Experimental Conditions --
sensitivity.pop <- rbind(update_data2(0,"rich", "Rich (0% CV)", "EM & UM"),
update_data2(5,"rich", "Rich (5% CV)", "EM & UM"),
update_data2(10,"rich", "Rich (10% CV)", "EM & UM"),
update_data2(20,"rich", "Rich (20% CV)", "EM & UM"),
update_data2(0,"sparse3pt", "Sparse (3pt, 0% CV)", "EM & UM"),
update_data2(5,"sparse3pt", "Sparse (3pt, 5% CV)", "EM & UM"),
update_data2(10,"sparse3pt", "Sparse (3pt, 10% CV)", "EM & UM"),
update_data2(20,"sparse3pt", "Sparse (3pt, 20% CV)", "EM & UM"),
update_data2(0,"sparse4pt", "Sparse (4pt, 0% CV)", "EM & UM"),
update_data2(5,"sparse4pt", "Sparse (4pt, 5% CV)", "EM & UM"),
update_data2(10,"sparse4pt", "Sparse (4pt, 10% CV)", "EM & UM"),
update_data2(20,"sparse4pt", "Sparse (4pt, 20% CV)", "EM & UM"),
update_data2(0,"PM", "Rich (0% CV)", "IM & PM"),
update_data2(5,"PM", "Rich (5% CV)", "IM & PM"),
update_data2(10,"PM", "Rich (10% CV)", "IM & PM"),
update_data2(20,"PM", "Rich (20% CV)", "IM & PM"),
update_data2(0,"PM_3pt", "Sparse (3pt, 0% CV)", "IM & PM"),
update_data2(5,"PM_3pt", "Sparse (3pt, 5% CV)", "IM & PM"),
update_data2(10,"PM_3pt", "Sparse (3pt, 10% CV)", "IM & PM"),
update_data2(20,"PM_3pt", "Sparse (3pt, 20% CV)", "IM & PM"),
update_data2(0,"PM_4pt", "Sparse (4pt, 0% CV)", "IM & PM"),
update_data2(5,"PM_4pt", "Sparse (4pt, 5% CV)", "IM & PM"),
update_data2(10,"PM_4pt", "Sparse (4pt, 10% CV)", "IM & PM"),
update_data2(20,"PM_4pt", "Sparse (4pt, 20% CV)", "IM & PM"))
sensitivity.pop$Condition <- factor(sensitivity.pop$Condition,
levels = unique(sensitivity.pop$Condition))Publication Figures
CYP2D6 Mediated 4-OH Atomoxetine Metabolism
4-OH Atomoxetine Formation (EM)
sensitivity_mid_end.plot <- ggplot(data = sensitivity.data %>%
filter(S <=100,
!Diplotype %in% c("1/2x2", "4/41", "4/5")),
aes(x = S, y = V, group = Diplotype))+
geom_line(data = sensitivity.pop %>%
filter(S<=100, !Diplotype %in% c("1/2x2", "4/41", "4/5")),
color = "grey",
size = 0.5,
aes(x = S, y = V, group = ID))+
geom_point(data = sensitivity.pop %>%
filter(S<=100, !Diplotype %in% c("1/2x2", "4/41", "4/5")),
aes(x = S, y = V, group = ID, color = Diplotype),
alpha = 0.4,
size = 2)+
geom_line(aes(color = Diplotype), size = 1)+
facet_wrap(~Condition, ncol = 4, nrow = 3, scales = "free")+
theme_bw()+
scale_color_viridis_d()+
xlab("Atomoxetine (uM)")+
ylab("4-OH-Atomoxetine Formation\n(pmol/min/mg protein)")+
ggeasy::easy_move_legend(to = "top")
sensitivity_mid_end.plot+
ggeasy::easy_add_legend_title("Genotype")
In-silico NLME predicted Michaelis-Menten fits of 4-hydroxy-atomoxetine formation using human liver microsomes by CYP2D6 extensive metabolizers (1/1, 1/2, 1/4, 2/2, 2/3, and 2/4) across all experimental conditions (sparse 3-4pt, 0-20% CV).
4-OH Atomoxetine Formation (UM)
# Rapid Metabolizers -------
sensitivity_high_end.plot <- ggplot(data = sensitivity.data %>%
filter(Diplotype == "1/2x2"),
aes(x = S, y = V, group = Diplotype))+
geom_line(data = sensitivity.pop %>%
filter(Diplotype == "1/2x2"),
color = "grey",
aes(x = S, y = V, group = ID),
size = 1)+
geom_point(data = sensitivity.pop %>%
filter(Diplotype == "1/2x2"),
aes(x = S, y = V, group = ID,color = Diplotype),
alpha = 0.4,
size = 2.5)+
geom_line(aes(color = Diplotype), size = 1.65)+
facet_wrap(Group~Condition, ncol = 4, nrow = 3, scales = "free")+
sjPlot::theme_sjplot()+
scale_color_viridis_d(option = "B", begin = 0)+
xlab("Atomoxetine (uM)")+
ylab("4-OH-Atomoxetine Formation\n(pmol/min/mg protein)")+
ggeasy::easy_move_legend(to = "top")
sensitivity_high_end.plot +
ggeasy::easy_add_legend_title("Genotype")
In-silico NLME predicted Michaelis-Menten fits of 4-hydroxy-atomoxetine formation using human liver microsomes by CYP2D6 ultra-rapid metabolizers (1/2X2) across all experimental conditions (sparse 3-4pt, 0-20% CV).
4-OH Atomoxetine Formation (IM/PM)
# Poor Metabolizers ----
sensitivity_low_end.plot <- ggplot(data = sensitivity.data %>%
filter(Diplotype %in% c("4/41", "4/5")),
aes(x = S, y = V, group = Diplotype))+
geom_line(data = sensitivity.pop %>%
filter(Diplotype %in% c("4/41", "4/5")),
color = "grey",
size = 1,
aes(x = S, y = V, group = ID))+
geom_line(aes(color = Diplotype), size = 1.65)+
geom_point(data = sensitivity.pop %>%
filter(Diplotype %in% c("4/41", "4/5")),
aes(x = S, y = V, group = ID,color = Diplotype),
alpha = 0.4,
size = 2.5)+
facet_wrap(Group~Condition, ncol = 4, nrow = 3, scales = "free")+
sjPlot::theme_sjplot()+
scale_color_viridis_d(option = "A", begin = 0.2, end = 0.5)+
# scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
# labels = trans_format("log10", math_format(10^.x)))+
xlab("Atomoxetine (uM)")+
ylab("4-OH-Atomoxetine Formation\n(pmol/min/mg protein)")+
ggeasy::easy_move_legend(to = "top")
sensitivity_low_end.plot+
ggeasy::easy_add_legend_title("Genotype")
In-silico NLME predicted Michaelis-Menten fits of 4-hydroxy-atomoxetine formation using human liver microsomes by CYP2D6 intermediat and poor metabolizers (4/41, and 4/5) across all experimental conditions (sparse 3-4pt, 0-20% CV).
CYP2D6 Intrinsic Clearance Predicted vs Actual
# Data frame with intrinsic clearance estimates for all conditions ----
sensitivity.CL <- sensitivity.data %>%
filter(S == 0) %>%
mutate(CLint = round(Vmax/Km,digits = 2)) %>%
select(-S, -V) %>%
left_join(diplotype.data %>%
mutate(CLint_actual = round(Vmax/Km, digit = 2)) %>%
select(Diplotype, CLint_actual),by = "Diplotype") %>%
mutate(CLint_ratio = CLint/CLint_actual)%>%
group_by(Diplotype) %>%
mutate(standard_resid = (CLint-CLint_actual)/sqrt(CLint_actual),
Obs = seq_along(standard_resid),
`Simulated Error` = if_else(Condition %in% c("Rich (0% CV)",
"Sparse (3pt, 0% CV)",
"Sparse (4pt, 0% CV)"),
"Control\n(0% CV)",
"Variable Error\n(5-20% CV)")) %>%
ungroup()
CL.plot <- ggplot(data = sensitivity.CL, aes(x = CLint_actual, y = CLint))+
geom_point(size = 3.5, alpha = 0.65, aes(fill = Condition), shape = 21)+
geom_abline(slope = 1, intercept = 0, color = "red")+
geom_abline(slope = 1, intercept = 0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.2, color = "red", linetype = "dashed")+
geom_abline(slope = 1, intercept = 0.15, color = "blue", linetype = "dashed")+
geom_abline(slope = 1, intercept = -0.15, color = "blue", linetype = "dashed")+
scale_y_log10()+
scale_x_log10()+
# facet_wrap(.~Condition, ncol = 4, scales = "free")+
xlab("Intrinsic Clearance (Actual)")+
ylab("Intrinsic Clearance (Predicted)")+
theme_bw(base_size = 14)+
ggeasy::easy_move_legend("right")+
scale_fill_viridis_d()+
# coord_cartesian(clip = "off")+
facet_wrap(~Group, ncol = 1, scales = "free")CLint Predicted vs Actual
CL.plot 
CYP2D6 Intrinsic Clearance Residual (1/2)
CL_plot2.1 <- ggplot(data = sensitivity.CL,
aes(y = standard_resid, x = Diplotype))+
geom_hline(yintercept = 0, size = 1, alpha = 0.5)+
geom_line(aes(color = Condition, group = Condition))+
geom_point(size = 3.5, shape = 21, alpha = 0.5, aes(fill = Condition))+
theme_bw(base_size = 14)+
facet_grid(`Simulated Error`~., switch = "both")+
scale_y_continuous(limits = c(-3,3))+
xlab("CYP2D6 Genotype")+
ylab(expression(atop("Standardized Residuals", "("*CL[intrinsic]*")")))
CL_plot2.1
CYP2D6 Intrinsic Clearance Residual (2/2)
CL_plot2.2 <- ggplot(data = sensitivity.CL,
aes(y = fct_reorder(Condition, standard_resid, .fun = median),
x = standard_resid))+
stat_boxplot(geom = "errorbar",
width = 0.5, size = 0.7) +
geom_boxplot(size = 0.7)+
geom_vline(xintercept = 0,
color = "blue", size = 1)+
geom_point(size = 2.5, shape = 21, aes(fill = Diplotype))+
# stat_summary(geom = "point",
# fun = mean,
# color = "red", size = 3)+
scale_x_continuous(limits = c(-3,3))+
theme_bw(base_size = 14)+
xlab(expression(atop("Standardized Residuals", "("*CL[intrinsic]*")")))+
ylab("Experimental Condition")+
geom_vline(xintercept = c(-2,2),
linetype = "dashed",
color = "darkgrey",
size = 1)+
ggeasy::easy_add_legend_title("CYP2D6\nGenotype")
CL_plot2.2
Diagnostic Plots
Residual Plots
Rich
plot(rich_covar.nlme, main = "Rich (0% CV)")
plot(rich_CV5_covar.nlme, main = "Rich (5% CV)")
plot(rich_CV10_covar.nlme, main = "Rich (10% CV)")
plot(rich_CV20_covar.nlme, main = "Rich (20% CV)")
Sparse 3pt
plot(sparse3pt_covar.nlme, main = "Sparse (3pt, 0% CV)")
plot(sparse3pt_CV5_covar.nlme, main = "Sparse (3pt, 5% CV)")
plot(sparse3pt_CV10_covar.nlme, main = "Sparse (3pt, 10% CV)")
plot(sparse3pt_CV20_covar.nlme, main = "Sparse (3pt, 20% CV)")
Sparse 4pt
plot(sparse4pt_covar.nlme, main = "Sparse (4pt, 0% CV)")
plot(sparse4pt_CV5_covar.nlme, main = "Sparse (4pt, 5% CV)")
plot(sparse4pt_CV10_covar.nlme, main = "Sparse (4pt, 10% CV)")
plot(sparse4pt_CV20_covar.nlme, main = "Sparse (4pt, 20% CV)")
PM Rich
plot(PM_covar.nlme, main = "PM Rich (0% CV)")
plot(PM_CV5_covar.nlme, main = "PM Rich (5% CV)")
plot(PM_CV10_covar.nlme, main = "PM Rich (10% CV)")
plot(PM_CV20_covar.nlme, main = "PM Rich (20% CV)")
PM Sparse 3pt
plot(PM_3pt_covar.nlme, main = "PM Sparse (3pt, 0% CV)")
plot(PM_3pt_CV5_covar.nlme, main = "PM Sparse (3pt, 5% CV)")
plot(PM_3pt_CV10_covar.nlme, main = "PM Sparse (3pt, 10% CV)")
plot(PM_3pt_CV20_covar.nlme, main = "PM Spase (3pt, 20% CV)")
PM Sparse 4pt
plot(PM_4pt_covar.nlme, main = "PM Sparse (4pt, 0% CV)")
plot(PM_4pt_CV5_covar.nlme, main = "PM Sparse (4pt, 5% CV)")
plot(PM_4pt_CV10_covar.nlme, main = "PM Sparse (4pt, 10% CV)")
plot(PM_4pt_CV20_covar.nlme, main = "PM Sparse (4pt, 20% CV)")
Predicted vs Observed Plots by CYP2D6 Genotype (UM/EM)
Rich Designs
plot(rich_covar.nlme,log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (0% CV)",
xlab = "log(Fitted values)")
plot(rich_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (5% CV)",
xlab = "log(Fitted values)")
plot(rich_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (10% CV)",
xlab = "log(Fitted values)")
plot(rich_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Rich (20% CV)",
xlab = "log(Fitted values)")
Sparse Designs (3pt)
plot(sparse3pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 0% CV)",
xlab = "log(Fitted values)")
plot(sparse3pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 5% CV)",
xlab = "log(Fitted values)")
plot(sparse3pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 10% CV)",
xlab = "log(Fitted values)")
plot(sparse3pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (3pt, 20% CV)",
xlab = "log(Fitted values)")
Sparse Designs (4pt)
plot(sparse4pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 0% CV )",
xlab = "log(Fitted values)")
plot(sparse4pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 5% CV)",
xlab = "log(Fitted values)")
plot(sparse4pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 10% CV)",
xlab = "log(Fitted values)")
plot(sparse4pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), main = "Sparse (4pt, 20% CV)",
xlab = "log(Fitted values)")
Predicted vs Observed Plots by CYP2D6 Genotype (IM/PM)
Rich Designs
plot(PM_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (0% CV)")
plot(PM_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (5% CV)")
plot(PM_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (10% CV)")
plot(PM_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Rich (20% CV)")
Sparse Designs (3pt)
plot(PM_3pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt, 0% CV)")
plot(PM_3pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt, 5% CV)")
plot(PM_3pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt 10% CV)")
plot(PM_3pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (3pt, 20% CV)")
Sparse Designs (4pt)
plot(PM_4pt_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 0% CV)")
plot(PM_4pt_CV5_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 5% CV)")
plot(PM_4pt_CV10_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 10% CV)")
plot(PM_4pt_CV20_covar.nlme, log(V)~log(fitted(.))|Diplotype,
abline = c(0,1), xlab = "log(Fitted values)",
main = "PM Sparse (4pt, 20% CV)")
Bootstrap Analysis
Virtual Population Function
An updated version of a previous function
population.sim() was generated called
population.sim2(). The updated function provides greater
flexibility substrate concentration range below 0.50 uM.
# ========== POPULATION SIMULATION FUNCTION =================#
## Updated to include low end point of 0.10 uM
population.sim2 <- function(
# Pre-define number of subjects per diplotype group
n = 10,
`CV%_D` = 25,
`CV%_V` = 0,
# Pre-define Michaelis Menten Kinetic Settings
Start = 0.5,
End = 2000,
By = 0.5,
# Turn off or on use of true population frequency
## If set to true "n" will equal total population number
pop_freq = FALSE,
# Pre-define reference data, and set seed
data = diplotype.data,
seed = 23457){
# ======== Create Data Containers ========
datasets <- list()
pop.data <- tibble()
# ======== Create Datasets by Diplotype ========
for (i in seq_along(data[[1]])){
# Setting Population Specific Frequency
n_pop <- if_else(pop_freq == FALSE, n,
ifelse(pop_freq == TRUE,
round(n*data[[i, c("Freq")]]), 0))
# Specify Diplotype for current loop
Diplotype = rep(data[[i,1]], n_pop)
# Random assignment of Vmax values N(0,sigma)
set.seed(seed)
Vmax_eta <- rnorm(n_pop, sd = `CV%_D`)
Vmax = data[[i,2]] + (Vmax_eta/100)*data[[i,2]]
# Random assignment of Km values N(0,sigma)
set.seed(seed)
Km_eta <- rnorm(n_pop, sd = `CV%_D`)
Km = data[[i,3]] + (Km_eta/100)*data[[i,3]]
# Store each diplotype group in a list container
temp <- tibble(Diplotype, Vmax, Km)
datasets[[i]] <- temp
}
# ======== Combine Diplotype Datasets ========
for (i in seq_along(data[[1]])){
pop.data <- rbind(pop.data, datasets[[i]])
}
# ======== Simulate Michaelis Menten ========
set.seed(seed)
pop.data <- pop.data %>%
mutate(ID = factor(seq_along(Diplotype),
levels = seq_along(Diplotype))) %>%
relocate(ID) %>%
expand_grid(S = c(0.1, seq(Start, End, By))) %>%
mutate(V = round((Vmax*S)/(Km+S), digits = 2))%>%
group_by(ID) %>%
# Additon of Residual Error
mutate(resid_pct = round(rnorm(n = length(S),sd = `CV%_V`),
digits = 3),
V_adj = V+(V*resid_pct/100)) %>%
ungroup()
return(pop.data)
}Generating Virtual Population (with Residual Error)
Michaelis-Menten kinetics were simulated for each subject in the virtual population across the 4 experimental designs (0%, 5%, 10%, and 20% CV).
# ========== CREATE BOOTSTRAP VIRTUAL POPULATIONS ============#
# Virtual Populations
CV0_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 0, End = 2000)
CV5_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 5, End = 2000)
CV10_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 10, End = 2000)
CV20_virtual.pop <- population.sim2(n = 1000, `CV%_V` = 20, End = 2000)
# Sample Population Containers
CV0_Sample_populations = list()
CV5_Sample_populations = list()
CV10_Sample_populations = list()
CV20_Sample_populations = list()Bootstrap Sampling of Virtual Population
Using a for() loop, a total of 10 bootstrapped
populations were generated for each experiental design (0-20% CV) by
randomly sampling (with replacement) 1% of the corresponding virtual
population (10 subject per genotype group, n = 90).
## Bootstrap for-loop
B = 10 ## number of bootstraps
for (i in 1:B) {
# Bootstrap for 0% CV Population ---------
CV0_Sample_populations[[i]] <- CV0_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
# Bootstrap for 5% CV Population ---------
CV5_Sample_populations[[i]] <- CV5_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
# Bootstrap for 10% CV Population ---------
CV10_Sample_populations[[i]] <- CV10_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
# Bootstrap for 20% CV Population ---------
CV20_Sample_populations[[i]] <- CV20_virtual.pop %>%
group_by(Diplotype) %>%
filter(ID %in% sample(ID, 10, replace = T)) %>%
ungroup()
}
# Combining Population Data -----------
Bootstrap.populations <- list(CV0_Sample_populations,
CV5_Sample_populations,
CV10_Sample_populations,
CV20_Sample_populations)
names(Bootstrap.populations)<- c("CV0_Pops",
"CV5_Pops",
"CV10_Pops",
"CV20_Pops")
# saveRDS(Bootstrap.populations, "R Output/Bootsrap Populations Data.rds")Sampling Function
Strategic sampling designs were implemented using a custom function
update_pop(), which functions similarly to the
update_data() function mentioned previously, however, is
designed to work more efficiently the Bootstrap.populations
data which is a list of dataframes). The inputs for the
update_pop() function include the following:
pop_data - A dataframe or indexed dataframe from a list of dataframes.
type - Specifies the strategic sampling design to be implemented. Options include: “rich”, “sparse3pt”, “sparse4pt”, “PM”, “PM_3pt”, “PM_4pt”.
#======== CREATE BOOTSTRAP EXPERIMENTAL DATASETS =========#
# Function to transform virtual population into sample population
update_pop <-function(pop_data, type){
# Sampling Information --------------
# Full Range Set
full_set <- c(0.1, 0.5, 1, 2.5, 5, 10, 25, 50, 100)
PM_range <- c(1, 5, 10, 15, 25, 30, 50, 60, 90, 100, 250, 500, 800, 1000, 1600, 2000)
# Strategic Sampling Sets ----
sparse3pt_set1 <- c(0.1, 2.5, 50)
sparse3pt_set2 <- c(1, 10, 100)
sparse4pt_set1 <- c(0.1, 2.5, 25, 50)
sparse4pt_set2 <- c(1, 10, 50, 100)
PM_3pt_set1 <- c(10,100,1000)
PM_3pt_set2 <- c(25,250,2000)
PM_4pt_set1 <- c(1,10,100,1000)
PM_4pt_set2 <- c(5,25,250,2000)
PM_range_set1 <- c(1,10,25,50,100,400,1000,2000)
PM_range_set2 <- c(5,15,30,60,90,200,800,1600)
ID_key <- pop_data %>%
select(Diplotype, ID) %>%
unique() %>%
group_by(Diplotype) %>%
mutate(n = 1:n()) %>%
ungroup()
PM_Diplotype <- ID_key %>%
select(Diplotype) %>%
filter(Diplotype %in% c("4/5","4/41")) %>%
unique()
EM_Diplotype <- ID_key %>%
select(Diplotype) %>%
filter(!Diplotype %in% c("4/5","4/41")) %>%
unique()
# Rich Sampling Simulation ----------------
Data <- pop_data %>%
# Filter S Based on Experimental Model Type
filter(S %in% switch(type,
"rich" = full_set,
"sparse3pt" = full_set,
"sparse4pt" = full_set,
"PM" = PM_range,
"PM_3pt" = PM_range,
"PM_4pt" = PM_range)) %>%
#Assign Sampling Key
left_join(ID_key, by = c("Diplotype", "ID")) %>%
# Strategic Sampling
mutate(Set = if_else(n %% 2 == 0, "Set 2", "Set 1")) %>%
filter(if_else(Set == "Set 1",
S %in% switch(type,
"rich" = full_set,
"sparse3pt" = sparse3pt_set1,
"sparse4pt" = sparse4pt_set1,
"PM" = PM_range_set1,
"PM_3pt" = PM_3pt_set1,
"PM_4pt" = PM_4pt_set1),
S %in% switch(type,
"rich" = full_set,
"sparse3pt" = sparse3pt_set2,
"sparse4pt" = sparse4pt_set2,
"PM" = PM_range_set2,
"PM_3pt" = PM_3pt_set2,
"PM_4pt" = PM_4pt_set2))) %>%
mutate(Diplotype = as_factor(Diplotype),
V = round(V_adj, digits = 3)) %>%
select(ID, Diplotype, S, V, Set) %>%
# Refine Inidvidual in dataset according to model
filter(Diplotype %in% switch(type,
"rich" = EM_Diplotype[[1]],
"sparse3pt" = EM_Diplotype[[1]],
"sparse4pt" = EM_Diplotype[[1]],
"PM" = PM_Diplotype[[1]],
"PM_3pt" = PM_Diplotype[[1]],
"PM_4pt" = PM_Diplotype[[1]]))
# Data Output ----------
Data
}In-Silico Experiments (Rich and Strategic Sampling)
Using a for() loop, each virtual population was subject
to all 3 sampling strategies (rich, sparse 3pt, and sparse 4pt); with
sampling range varying according to CYP2D6 genotype as described
previously (see Strategic Sampling Approach). Data was combined
as a list of a list of dataframes with the first level (n = 4) being %CV
specific data, the second level (n = 3) being sampling strategy with
each variable at this level containing data for 10 unique population (n
= 120 total datasets).
# Create Experimental Data sets Populations ---------
CV0_Boot.data <- list()
CV5_Boot.data <- list()
CV10_Boot.data <- list()
CV20_Boot.data <- list()
# For Loop to apply strategic sampling to all dataframes ----------
for (i in 1:B) {
CV0_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "rich")
CV0_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "sparse3pt")
CV0_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "sparse4pt")
CV0_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "PM")
CV0_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "PM_3pt")
CV0_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV0_Pops[[i]], type = "PM_4pt")
cat('CV0 Iteration', i, "of", B, "complete...\n")
CV5_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "rich")
CV5_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "sparse3pt")
CV5_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "sparse4pt")
CV5_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "PM")
CV5_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "PM_3pt")
CV5_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV5_Pops[[i]], type = "PM_4pt")
cat('CV5 Iteration', i, "of", B, "complete...\n")
CV10_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "rich")
CV10_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "sparse3pt")
CV10_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "sparse4pt")
CV10_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "PM")
CV10_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "PM_3pt")
CV10_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV10_Pops[[i]], type = "PM_4pt")
cat('CV10 Iteration', i, "of", B, "complete...\n")
CV20_Boot.data$rich[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "rich")
CV20_Boot.data$sparse3pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "sparse3pt")
CV20_Boot.data$sparse4pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "sparse4pt")
CV20_Boot.data$PM[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "PM")
CV20_Boot.data$PM_3pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "PM_3pt")
CV20_Boot.data$PM_4pt[[i]] <- update_pop(Bootstrap.populations$CV20_Pops[[i]], type = "PM_4pt")
cat('CV20 Iteration', i, "of", B, "complete...\n")
}
# Combine all boostrap data into a single data structure ---------
Complete_Bootstrap.data <- list(CV0_Boot.data,
CV5_Boot.data,
CV10_Boot.data,
CV20_Boot.data)
names(Complete_Bootstrap.data) <-c("CV0","CV5","CV10","CV20")
# Save the completed bootstrap data sets ------
# saveRDS(Complete_Bootstrap.data, "R Output/Complete Bootstrap Data.rds")NLME Fitting Function
fit_nlme() is a custom function used to simultaneously
fit non-linear least squares, and non-linear mixed effect models (with
and without covariates). The inputs for the function include the
following:
data - Data to be used for model fitting
type - Type of CYP2D6 metabolizer. Options include: “EM” for extensive metabolizer and ultra-rapid metabolizer data, and “PM” for intermediate and poor metabolizer data (default = “EM”).
which - Specifies the desired model output. Options include: “nls” = non-linear least squares, “base.model” for NLME with out covariates, “covar.model” for NLME with covariates (CYP2D6 genotype); default = “covar.model”.
# ==== Fitting Mixed Effect Model ================#
fit_nlme <- function(data, type = "EM", which = "covar.model"){
n <- if_else(type == "PM", 1, 6)
# Initial Screen ----
data0.grp <- groupedData(V~S|ID, data)
fit0.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = data0.grp)
# Extract Erroneous Subjects
error.vmax <- coef(fit0.nls) %>%
filter(is.na(Vmax)) %>%
row.names() %>%
as.numeric()
error.km <- coef(fit0.nls) %>%
filter(is.na(Km)) %>%
row.names() %>%
as.numeric()
error.list <- c(error.vmax, error.km)
# Clean Data set
data.grp <- data0.grp %>% filter(!ID %in% error.list)
# Base Model Fits
fit.nls <- nlsList(V~SSmicmen(S,Vmax,Km), data = data.grp)
base.nlme <- nlme(fit.nls, random = pdDiag(Vmax + Km ~ 1))
covar.nlme <- update(base.nlme, fixed = Vmax + Km ~ Diplotype,
start = c(fixef(base.nlme)[[1]], rep(0,n),
fixef(base.nlme)[[2]], rep(0,n)),
weights = varPower())
# Output ----------
switch (which,
"nls" = fit.nls,
"base.model" = base.nlme,
"covar.model" = covar.nlme)
}NLME Estimate Extraction Function
extract_nlme_est() is a custom function used to fit,
extract, and tidy NLME parameter estimates from experimental data. This
function is optimized to work within the context of a for()
loop, and contains the fit_nlme() function as a dependency.
The inputs for the function include the following:
exp.data - Experimental data that NLME estimates are desired from.
condition - Label used to tag experimental condition from which the experimental data belongs.
PM - TRUE or FALSE argument specifying where data belongs to IM/PM CYP2D6 genotype (TRUE), or UM/EM CYP2D6 genotype (FALSE); default = “FALSE”.
pop.num - Label specifying which bootstrapped population is being analyzed (default = 1).
extract_nlme_est <- function(exp.data, condition, PM = FALSE, pop.num = 1){
# Input ----------
model.type <- if_else(PM == TRUE, "PM", "EM")
nlme.model <- fit_nlme(exp.data, type = model.type)
int.data <- intervals(nlme.model, which = "fixed")
label <- if_else(PM == TRUE, "4/41", "1/1")
# Tidy Data -------
temp <- int.data$fixed %>%
as_tibble() %>%
mutate(Parameter = rownames(int.data$fixed)) %>%
separate(Parameter, remove = T, sep = "\\.", into = c("Parameter", "Diplotype"))%>%
mutate(Diplotype = recode(Diplotype, "(Intercept)" = paste0("Diplotype",label))) %>%
group_by(Parameter) %>%
mutate(est. = if_else(Diplotype != paste0("Diplotype",label),
est. + est.[Diplotype == paste0("Diplotype",label)],
est.),
across(where(is.numeric), round, 2)) %>%
ungroup() %>%
separate(Diplotype, remove = T, sep = "Diplotype", into = c(".", "Diplotype"))%>%
select(Parameter, Diplotype, est.)
colnames(temp) <- c("Parameter", "Diplotype", paste(condition))
# Output ----
temp %>%
pivot_longer( 3, names_to = "Condition", values_to = "Estimate") %>%
mutate(Population = pop.num)
}
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$sparse3pt[[3]],
condition = "Sparse (3pt, CV 5%)", PM = FALSE, pop.num = 3)Non-Linear Mixed Effect Modeling (Bootstrap Populations)
Using a for() loop, bootstrapped parameter estimates
were extracted from experimental data
(Complete_Bootstrap.data) using the
extract_nlme.est() function. Parameter estimate for all
experimental designs and conditions for each population were stored in
object Complete_Bootstrap.est as a list. A collapsed
version was saved as Complete Bootstrap table.RDS for data
analysis and user reproducibility purposes.
Complete_Bootstrap.est <- list()
for (i in 1:B) {
# CV 0% Condition -------------
CV0_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$rich[[i]],
condition = "Rich (CV 0%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$sparse3pt[[i]],
condition = "Sparse (3pt, CV 0%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$sparse4pt[[i]],
condition = "Sparse (4pt, CV 0%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$PM[[i]],
condition = "Rich (CV 0%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$PM_3pt[[i]],
condition = "Sparse (3pt, CV 0%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV0$PM_4pt[[i]],
condition = "Sparse (4pt, CV 0%)", PM = TRUE, pop.num = i))
cat('CV0% - Iteration', i, "of", B, "complete...\n")
# # # CV 5% Condition -------------
CV5_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$rich[[i]],
condition = "Rich (CV 5%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$sparse3pt[[i]],
condition = "Sparse (3pt, CV 5%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$sparse4pt[[i]],
condition = "Sparse (4pt, CV 5%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$PM[[i]],
condition = "Rich (CV 5%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$PM_3pt[[i]],
condition = "Sparse (3pt, CV 5%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV5$PM_4pt[[i]],
condition = "Sparse (4pt, CV 5%)", PM = TRUE, pop.num = i))
cat('CV5% - Iteration', i, "of", B, "complete...\n")
# # CV 10% Condition -------------
CV10_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$rich[[i]],
condition = "Rich (CV 10%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$sparse3pt[[i]],
condition = "Sparse (3pt, CV 10%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$sparse4pt[[i]],
condition = "Sparse (4pt, CV 10%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$PM[[i]],
condition = "Rich (CV 10%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$PM_3pt[[i]],
condition = "Sparse (3pt, CV 10%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV10$PM_4pt[[i]],
condition = "Sparse (4pt, CV 10%)", PM = TRUE, pop.num = i))
cat('CV10% - Iteration', i, "of", B, "complete...\n")
# CV 20% Condition -------------
# Problematic fails to converge, requires further attention ##
CV20_est.table <- rbind(
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$rich[[i]],
condition = "Rich (CV 20%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$sparse3pt[[i]],
condition = "Sparse (3pt, CV 20%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$sparse4pt[[i]],
condition = "Sparse (4pt, CV 20%)", pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$PM[[i]],
condition = "Rich (CV 20%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$PM_3pt[[i]],
condition = "Sparse (3pt, CV 20%)", PM = TRUE, pop.num = i),
extract_nlme_est(exp.data = Complete_Bootstrap.data$CV20$PM_4pt[[i]],
condition = "Sparse (4pt, CV 20%)", PM = TRUE, pop.num = i))
cat('CV20% - Iteration', i, "of", B, "complete...\n")
Complete_Bootstrap.est[[i]] <- rbind(CV0_est.table,
CV5_est.table,
CV10_est.table,
CV20_est.table)
# Progress Report -------
cat('==== ROUND', i, "of", B, "COMPLETE! ===\n")
}
Complete_Bootstrap.table <- bind_rows(Complete_Bootstrap.est)
# saveRDS(R Output/Complete Bootstrap Table.rds")Bootstrap Evaluation
Non-parametric bootstrap analysis based on a total of 120 generated
data sets (10 re-sampled populations per experimental condition) was
performed to evaluate the precision of the final model parameters (see
separate Bootstrap Analysis section above). Mean parameter
estimates and 95% confidence intervals for each condition were
summarized by genotype. Confidence intervals were calculated using the a
custom confidence interval function CI() detailed below
based on (Equation 6); where confidence interval (CI) is equal to the
mean value of the sample population plus or minus the t-score at
significance level(0.05) for N-1 degrees of freedom multiplied by the
standard error of the sample population mean. The inputs for the
CI() function are as followed:
x - A vector of values for which a confidence interval is to be calculated.
alpha - Alpha level for confidence interval (default = 0.05)
which - Specifies the desired output. Options: confidence interval rang = “ci”, lower interval estimate = “lower”, upper interval estimate = “upper”, mean estimate = “est”, all estimates (mean [lower - upper]) = “all”.
Confidence Interval (95%) Function
# Import Bootstrap Data ----------
Bootstrap.table <- read_rds("R Output/Complete Bootstrap Table.rds")
# function to calculate confidence interval of bootstrap estimates -----
CI <- function(x, alpha = 0.05, which = "ci"){
mean <- mean(x)
n <- length(x)
std_dev <- sd(x)
std_error <- std_dev/sqrt(n)
degrees_of_freedom = n - 1
t_score <- qt(p=alpha/2, df=degrees_of_freedom, lower.tail=F)
margin_error <- t_score * std_error
lower <- mean - margin_error
upper <- mean + margin_error
ci <- paste0("(", round(lower, 2),
" - ",round(upper, 2),")")
all <- paste0(round(mean, 2),
" (", round(lower, 2),
" - ",round(upper, 2),")")
# Output ------
switch (which,
"ci" = ci,
"lower" = lower,
"upper" = upper,
"est" = mean,
"all" = all)
}Bootstrap Analysis Summary Datatable
Summary.boot <- Bootstrap.table %>%
mutate(Index = 1:length(Population)) %>%
filter(Index != 315) %>% # un-physiologic scenario
group_by(Parameter, Diplotype, Condition) %>%
summarize(Est = round(mean(Estimate), 2),
sd = round(sd(Estimate), 2),
Lower = CI(Estimate, which = "lower"),
Upper = CI(Estimate, which = "upper"),
`CI (95%)` = CI(Estimate, which = "ci"),
all = CI(Estimate, which = "all"),
n = length(Estimate)) %>%
ungroup() %>%
left_join(actual.est, by = c("Diplotype", "Parameter")) %>%
mutate(Residual = Actual - Est)Bootstrap Estimates
Vmax Estimates (w/95% CI) Plot
# Vmax Plot --------------
ggplot(data = Summary.boot %>%
filter(Parameter == "Vmax"),
aes(x = Est, y = Condition))+
geom_vline(aes(xintercept = Actual),
color = "blue", size = 0.7, alpha = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.25)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*0.8)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.3)), alpha = 0)+
geom_vline(aes(xintercept = c(Actual*0.7)), alpha = 0)+
geom_errorbar(aes(xmin = Lower, xmax = Upper), width = 0.6)+
geom_pointrange(aes(xmin = Lower, xmax = Upper))+
facet_wrap(~Diplotype, scales = "free_x", ncol = 3, nrow = 3)+
theme_bw()+
ggtitle("Bootstrapped Confidence Intervals (95%) of Vmax")+
xlab(bquote(paste('V'['max']*" Estimate")))+
ylab("Experimental Condition")
Bootstrapped parameter estimates with 95% confidence intervals (CI) for final model across all experimental conditions (stratified by CYP2D6 genotype). Black markers with intervals represent the mean estimate and 95% confidence interval of 10 bootstrapped sample populations, blue line represents true estimate of Vmax for specified genotype group, dotted red lines represent 80%-125% lower-and-upper bounds of parameter estimate.
Km Estimates (w/95% CI) Plot
# Km Plot ---------
ggplot(data = Summary.boot %>%
filter(Parameter == "Km"),
aes(x = Est, y = Condition))+
geom_vline(aes(xintercept = Actual),
color = "blue", size = 0.7, alpha = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.25)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*0.8)), color = "red",
linetype = "dashed", size = 0.7)+
geom_vline(aes(xintercept = c(Actual*1.3)), alpha = 0)+
geom_vline(aes(xintercept = c(Actual*0.7)), alpha = 0)+
geom_errorbar(aes(xmin = Lower, xmax = Upper), width = 0.6)+
geom_pointrange(aes(xmin = Lower, xmax = Upper))+
facet_wrap(~Diplotype, scales = "free_x", ncol = 3, nrow = 3)+
theme_bw()+
ggtitle("Bootstrapped Confidence Intervals (95%) of Km")+
xlab(bquote(paste('K'['m']*" Estimate")))+
ylab("Experimental Condition")
Bootstrapped parameter estimates with 95% confidence intervals (CI) for final model across all experimental conditions (stratified by CYP2D6 genotype). Black markers with intervals represent the mean estimate and 95% confidence interval of 10 bootstrapped sample populations, blue line represents true estimate of Km for specified genotype group, dotted red lines represent 80%-125% lower-and-upper bounds of parameter estimate
R Session Information
# Session Information ------
sessionInfo()## R version 4.2.0 (2022-04-22 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kableExtra_1.3.4 cowplot_1.1.1 scales_1.2.0 ggpubr_0.4.0
## [5] ggeasy_0.1.3 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.9
## [9] purrr_0.3.4 readr_2.1.2 tidyr_1.2.0 tibble_3.1.8
## [13] ggplot2_3.3.6 tidyverse_1.3.1 nlme_3.1-157
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.2 lubridate_1.8.0 insight_0.17.1 webshot_0.5.4
## [5] httr_1.4.3 tools_4.2.0 backports_1.4.1 bslib_0.4.0
## [9] sjlabelled_1.2.0 utf8_1.2.2 R6_2.5.1 DBI_1.1.2
## [13] colorspace_2.0-3 withr_2.5.0 tidyselect_1.1.2 gridExtra_2.3
## [17] emmeans_1.7.4-1 compiler_4.2.0 performance_0.9.0 cli_3.3.0
## [21] rvest_1.0.2 xml2_1.3.3 labeling_0.4.2 bayestestR_0.12.1
## [25] sass_0.4.2 mvtnorm_1.1-3 systemfonts_1.0.4 digest_0.6.29
## [29] minqa_1.2.4 rmarkdown_2.14 svglite_2.1.0 pkgconfig_2.0.3
## [33] htmltools_0.5.2 lme4_1.1-29 dbplyr_2.2.0 fastmap_1.1.0
## [37] highr_0.9 rlang_1.0.4 readxl_1.4.0 rstudioapi_0.13
## [41] jquerylib_0.1.4 farver_2.1.1 generics_0.1.2 jsonlite_1.8.0
## [45] car_3.1-0 sjPlot_2.8.10 magrittr_2.0.3 Matrix_1.4-1
## [49] parameters_0.18.1 Rcpp_1.0.8.3 munsell_0.5.0 fansi_1.0.3
## [53] abind_1.4-5 lifecycle_1.0.1 stringi_1.7.6 yaml_2.3.5
## [57] carData_3.0-5 MASS_7.3-56 grid_4.2.0 sjmisc_2.8.9
## [61] crayon_1.5.1 lattice_0.20-45 splines_4.2.0 ggeffects_1.1.2
## [65] haven_2.5.0 sjstats_0.18.1 hms_1.1.1 knitr_1.39
## [69] pillar_1.8.0 boot_1.3-28 estimability_1.3 ggsignif_0.6.3
## [73] effectsize_0.7.0 codetools_0.2-18 reprex_2.0.1 glue_1.6.2
## [77] evaluate_0.15 modelr_0.1.8 nloptr_2.0.3 vctrs_0.4.1
## [81] tzdb_0.3.0 cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [85] datawizard_0.4.1 cachem_1.0.6 xfun_0.31 xtable_1.8-4
## [89] broom_0.8.0 coda_0.19-4 rstatix_0.7.0 viridisLite_0.4.0
## [93] ellipsis_0.3.2References
(1)
R
Core Team. R: A Language and
Environment for Statistical Computing; R Foundation for
Statistical Computing: Vienna, Austria, 2022.
(2)
Wilke, C. O. Cowplot: Streamlined Plot Theme and
Plot Annotations for Ggplot2; 2020.
(3)
Wickham, H.; François, R.; Henry, L.; Müller,
K. Dplyr: A
Grammar of Data Manipulation; 2022.
(4)
Carroll, J.; Schep, A.; Sidi, J. Ggeasy: Easy Access to
Ggplot2 Commands; 2021.
(5)
Wickham, H.; Chang, W.; Henry, L.; Pedersen, T.
L.; Takahashi, K.; Wilke, C.; Woo, K.; Yutani, H.; Dunnington, D. Ggplot2: Create
Elegant Data Visualisations Using the Grammar of Graphics;
2022.
(6)
Kassambara, A. Ggpubr: Ggplot2 Based
Publication Ready Plots; 2020.
(7)
(8)
Pinheiro, J.; Bates, D.; R Core Team. Nlme: Linear and
Nonlinear Mixed Effects Models; 2022.
(9)
Wickham, H.; Hester, J.; Bryan, J. Readr: Read Rectangular
Text Data; 2022.
(10)
Wickham, H.; Seidel, D. Scales: Scale Functions
for Visualization; 2022.
(11)
Müller, K.; Wickham, H. Tibble: Simple Data
Frames; 2022.
(12)
Wickham, H.; Girlich, M. Tidyr: Tidy Messy
Data; 2022.
(13)
Wickham, H. Tidyverse: Easily
Install and Load the Tidyverse; 2021.
(14)
Wickham, H. Ggplot2: Elegant Graphics for Data
Analysis; Springer-Verlag New York, 2016.
(15)
Pinheiro, J. C.; Bates, D. M. Mixed-Effects
Models in s and s-PLUS; Springer: New York, 2000. https://doi.org/10.1007/b98882.
(16)
Wickham, H.; Averick, M.; Bryan, J.; Chang, W.;
McGowan, L. D.; François, R.; Grolemund, G.; Hayes, A.; Henry, L.;
Hester, J.; Kuhn, M.; Pedersen, T. L.; Miller, E.; Bache, S. M.; Müller,
K.; Ooms, J.; Robinson, D.; Seidel, D. P.; Spinu, V.; Takahashi, K.;
Vaughan, D.; Wilke, C.; Woo, K.; Yutani, H. Welcome to the tidyverse. Journal of Open Source
Software 2019, 4 (43), 1686. https://doi.org/10.21105/joss.01686.